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Abstract 

Quantum hydrodynamics in superfluid helium and atomic Bose-Einstein con- 
densates (BECs) has been recently one of the most important topics in low 
temperature physics. In these systems, a macroscopic wave function (or- 
der parameter) appears because of Bose-Einstein condensation, which cre- 
ates quantized vortices. Turbulence consisting of quantized vortices is called 
quantum turbulence (QT). The study of quantized vortices and QT has in- 
creased in intensity for two reasons. The first is that recent studies of QT 
are considerably advanced over older studies, which were chiefly limited to 
thermal counterflow in 4 He, which has no analogue with classical traditional 
turbulence, whereas new studies on QT are focused on a comparison between 
QT and classical turbulence. The second reason is the realization of atomic 
BECs in 1995, for which modern optical techniques enable the direct control 
and visualization of the condensate and can even change the interaction; such 
direct control is impossible in other quantum condensates like superfluid he- 
lium and superconductors. Our group has made many important theoretical 
and numerical contributions to the field of quantum hydrodynamics of both 
superfluid helium and atomic BECs. In this article, we review some of the 
important topics in detail. The topics of quantum hydrodynamics are di- 
verse, so we have not attempted to cover all these topics in this article. We 
also ensure that the scope of this article does not overlap with our recent re- 
view article ( 1arXiv: 1004.5458|). "Quantized vortices in superfluid helium and 
atomic Bose-Einstein condensates" , and other review articles. 
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1. Introduction 

Quantum hydrodynamics (QHD) refers to the hydrodynamics of quan- 
tum condensed fluids, such as superfluid helium and atomic Bose-Einstein 
condensates (BECs), which are subject to quantum restrictions. As a re- 
sult of Bose-Einstein condensation, a system exhibits a macroscopic wave 
function ty(r,t) = \^(r, t)\e % ^ r ^ as an order parameter. The superfluid ve- 
locity field of a quantum condensed fluid is given by v s = (hfM)'V(f), with 
boson mass M, representing the potential flow. Since the macroscopic wave 
function should be single-valued for the space coordinate r, the circulation 
§ v ■ d£ for an arbitrary closed loop in the fluid is quantized by the quantum 
k = h/M. A vortex with such quantized circulation is called a quantized 
vortex. Any rotational motion of a superfluid is sustained only by quantized 
vortices. 

The appearance of quantized vortices distinguishes the hydrodynamics of 
quantum condensed fluids from that of classical fluids. A quantized vortex 
is a stable topological defect characteristic of a BEC and is different from 
a vortex in a classical viscous fluid. Firstly, the circulation is quantized, 
which contrasts with a classical vortex that can have any value of circulation. 
Secondly, a quantized vortex is a vortex of inviscid superflow and thus it 
cannot decay by the viscous diffusion of vorticity that occurs in a classical 
fluid. Thirdly, the core of a quantized vortex is very thin, on the order of the 
coherence length, which is only a few angstroms in superfluid 4 He and sub 
/jm even in atomic BECs. Since the vortex core is very thin and does not 
decay by diffusion, it is always possible to identify the position of a quantized 
vortex in the fluid. 

The turbulence of a superfluid velocity field is called superfluid turbu- 
lence or quantum turbulence (QT) [l], 0]. Since any rotational motion of a 
superfluid is sustained by quantized vortices, QT usually takes the form of a 
disordered tangle of quantized vortices. QT is currently the most important 
subject of research of QHD in the field of low temperature physics. The tur- 
bulence of classical fluids, called classical turbulence (CT), has been studied 
intensively in a number of fields, but it is still not yet well understood [3[. 
This is chiefly because turbulence is a complicated dynamical phenomenon 
with strong nonlinearity. Vortices may be the key to understanding turbu- 
lence. For example, Leonardo da Vinci observed the turbulent flow of water 
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and drew some sketches showing that turbulence had a structure comprised 
of vortices of different sizes. However, vortices are not well-defined for a clas- 
sical viscous fluid. They are unstable and appear and disappear repeatedly. 
The circulation is not conserved and not identical for each vortex. Comparing 
QT and CT reveals definite differences, which demonstrates the importance 
of studying QT. QT consists of a tangle of quantized vortices that have the 
same conserved circulation. Thus, QT can be an easier system to study than 
CT and present a much simpler model of turbulence than CT. 

QHD was first studied in superfluid 4 He and more recently in atomic 
BECs. Liquid 4 He enters a superfluid state below the A point (2.17 K) 



with Bose-Einstein condensation of the 4 He atoms |4|. The characteristic 
phenomena of superfluidity were discovered experimentally in the 1930s by 
Kapitza et al. The hydrodynamics of superfluid helium is well described by 
the two-fluid model, for which the system consists of an inviscid superfluid 
(density p s ) and a viscous normal fluid (density p n ) with two independent 
velocity fields v s and v n . The mixing ratio of the two fluids depends on 
temperature. As the temperature is reduced below the A point, the ratio of 
the superfluid component increases, and the entire fluid becomes a superfluid 
below approximately 1 K. Early experimental studies on superfluid turbu- 
lence focused primarily on thermal counterflow, in which the normal fluid 
and superfluid flow in opposite directions. The flow is driven by an injected 
heat current and it was found that the superflow becomes dissipative when 
the relative velocity between the two fluids exceeds a critical value. This was 
nothing but the appearance of QT. The interaction between the cores of tan- 
gled vortices and the normal fluid causes the dissipation. Considerable effort 
has been devoted to the study of thermal counterflow in superfluid 4 He. How- 
ever, since counterflow turbulence has no classical analog, the relationship 
between QT and CT has not been satisfactorily studied. 

Research into QHD has tended toward new directions since the mid 1990s. 
One new direction is in the field of low temperature physics, studying super- 
fluid helium. This field of study started with the attempt to understand the 
relationship between QT and CT j^l, @, 0]. The energy spectrum of fully 
developed CT is known to obey the Kolmogorov law in the inertial range. 
Recent experimental and numerical studies support a Kolmogorov spectrum 
in QT. Following these studies, QT research on superfluid helium has pro- 
gressed to important topics such as the dissipation process at very low tem- 
peratures, QT created by vibrating structures, and the visualization of QT 
0|. Another new direction is the realization of Bose-Einstein condensa- 
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tion in trapped atomic gases, first performed in 1995, which has stimulated 
intense experimental and theoretical activity js), [^j. As proof of the existence 
of superfluidity, quantized vortices have been created and observed in atomic 
BECs, and considerable effort has been devoted to a number of fascinating 



problems related to this [10|, [llj • Atomic BECs have several advantages over 
superfluid helium. The most important is that modern optical techniques 
enable the direct control of condensates and the visualization of quantized 
vortices. 

This article reviews the recent theoretical and numerical contributions 
of the Osaka City University group on QHD. Section 2 describes the basics 
of QHD. Section 3 describes QT using the vortex filament model and the 
Gross-Pitaevskii model. Section 4 describes hydrodynamic instability in two- 
component BECs. Section 5 is devoted to conclusions. 

2. Basics of quantum hydrodynamics 

This section reviews briefly the background of low temperature physics 
necessary for understanding this article. 

2.1. Bose-Einstein condensation 

Quantum condensation appearing in quantum fluids is caused chiefly by 
Bose-Einstein condensation. Quantum mechanics is often thought to give the 
physical laws at microscopic scales, but quantum mechanics appears even at 
macroscopic scales through BECs. 

The essence of quantum mechanics is the duality of the particle picture 
and wave picture. Let us consider an ideal atomic gas. At relatively high tem- 
peratures, the statistics of the atoms obey the classical Maxwell-Boltzmann 
distribution and each atom behaves like a particle. As the temperature is 
reduced, however, the thermal de Broglie wavelength is increased to become 
comparable to the mean distance between atoms. Then each atom begins to 
behave like a wave and the statistics changes to the quantum Fermi-Dirac or 
Bose-Einstein distribution depending on whether the atom is a Fermion or 
a Boson. If the atoms are Bosons and the system is cooled below a critical 
temperature T BEC , Bose-Einstein condensation occurs in which the atoms 
occupy the same single-particle ground state. The critical temperature is 
given by 

ft 2 " 2 ' 3 
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where the relevant quantities are the particle mass M, the number density n, 
the Planck constant h = 2nh, and the Boltzmann constant ks- Then, matter 
waves of atoms become coherent to form a macroscopic wave function (the 
order parameter) ^(r, t) = \^(r, t)\e ld ^ r '^ extending over the whole volume 
of the system, and the assemblage of these atoms is called a Bose-Einstein 
condensate (BEC). 

2.2. Liquid helium, superfluidity and the two- fluid model 

Independent of studies of quantum statistical mechanics, the field of low 
temperature physics has developed since the beginning of 20th century. Low 
temperature physics is generally believed to start with the first liquefaction 
of 4 He at 4.2 K by Onnes in 1908. Subsequently, Onnes observed supercon- 
ductivity in mercury in 1911. Onnes noticed an anomaly in the heat capacity 
of liquid helium at the A point T\ = 2.17 K. In 1938 Kapitza et al. observed 
that liquid 4 He becomes inviscid below the A point and called this striking 



phenomenon superfluidity [12j, [l3| . The superfluid transition of liquid 4 He at 
2.17 K is called the A transition. 

In order to explain the various hydrodynamic phenomena of superfluid- 
ity, Tisza [IH and Landau [l5| introduced the two-fluid model. According 
to the two-fluid model, the system consists of an inviscid superfluid (density 
p s ) and a viscous normal fluid (density p n ) with two independent velocity 
fields v s and v n . The mixing ratio of the two fluids depends on temper- 
ature. As the temperature is reduced below the A point, the ratio of the 
superfluid component increases and the entire fluid becomes a superfluid be- 
low approximately 1 K. While the two-fluid model successfully explained the 
phenomena of superfluidity, it was discovered in the 1940s that superfluidity 



breaks down when a superfluid flows fast [16( and this phenomenon could 
not be explained through the two-fluid model. This was later found to be 
caused by turbulence of the superfluid component due to random motion of 
quantized vortices. 

2.3. Bose-Einstein condensation, macroscopic wave function, and quantized 
vortices 

The A transition is closely related to the Bose-Einstein condensation of 
4 He atoms. London proposed theoretically in 1938 that the A transition of 
liquid 4 He is caused by Bose-Einstein condensation of 4 He atoms [l7|. When 
Tbec °f Eq. ((Tj) is evaluated for the mass and density appropriate to liquid 
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4 He at saturated vapor pressure, Tbec of approximately 3.13 K is obtained, 
which is close to T\ = 2.17 K. 

A Bose-condensed system exhibits the macroscopic wave function ty(r,t) = 
\^(r, t)\e l ^ r '^ as an order parameter. The superfluid velocity field is given 
by v s = (hfM)'V(f), representing the potential flow. Since the macroscopic 
wave function should be single-valued for the space coordinate r, the circu- 
lation § v ■ di for an arbitrary closed loop in the fluid is quantized by the 
quantum k = h/M. A vortex with quantized circulation is called a quantized 
vortex. Any rotational motion of a superfluid is sustained only by quantized 
vortices. 



3. Quantum turbulence 



3.1. Research history 

This subsection describes briefly the research history of quantized vortices 
and QT in superfluid 4 He. 

The idea of quantized circulation was first proposed by Onsager for a 



series of annular rings in a rotating superfluid |18|. Feynman considered 



that a vortex in a superfluid can take the form of a vortex filament with 
quantized circulation k and a core of atomic dimensions Early exper- 

imental studies on superfluid hydrodynamics focused primarily on thermal 
counterflow. The flow is driven by an injected heat current, and the normal 
fluid and superfluid flow in opposite directions. The superflow was found 
to become dissipative when the relative velocity between the two fluids ex- 
ceeds a critical value [l6j]. Gorter and Mellink attributed the dissipation to 
mutual friction between two fluids and considered the possibility of super- 
fluid turbulence. Feynman proposed a turbulent superfluid state consisting 



of a tangle of quantized vortices [191 ] . Hall and Vinen performed experiments 
of second sound attenuation in rotating 4 He, where second sound refers to 
the entropy wave in which superfluid and normal fluid oscillate oppositely, 
and its propagation and attenuation give information on the vortex density 
in the fluid. They found that mutual friction arises from the interaction 



between the normal fluid and quantized vortices [20|, |21|. Vinen confirmed 
Feynman's findings experimentally by showing that the dissipation in ther- 
mal counterflow arises from mutual friction between vortices and the normal 



flow [22j, |23j, |24], |25| . Vinen also succeeded in observing quantized circulation 
using vibrating wires in rotating superfluid 4 He 26]. Subsequently, many 



experimental studies have examined superfluid turbulence (ST) in thermal 
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counterflow systems and have revealed a variety of physical phenomena [27 . 
Since the dynamics of quantized vortices are nonlinear and non-local, it has 
not been easy to quantitatively understand these observations on the basis 
of vortex dynamics. Schwarz clarified the picture of ST based on tangled 
vortices by numerical simulation of the quantized vortex filament model in 



the thermal counterflow [281. 129j. However, since the thermal counterflow has 
no analogy in conventional fluid dynamics, this study was not helpful in clar- 
ifying the relationship between ST and classical turbulence (CT). ST is often 
called quantum turbulence (QT), which emphasizes the quantum effects. 

QHD, including QT, is reduced to the motion of quantized vortices. 
Hence, understanding the dynamics of quantized vortices is a key issue in 
QHD. Two formulations are generally available for studying the dynamics of 
quantized vortices. One is the vortex filament model and the other is the 
Gross-Pitaevskii (GP) model. This section describes the results first by the 
vortex filament model and then by the GP model. 

3.2. Vortex filament model 

As described in Sec. 1., a quantized vortex has quantized circulation. The 
vortex core is extremely thin, usually much smaller than other characteristic 
scales of vortex motion. These properties allow a quantized vortex to be 
represented as a vortex filament. In classical fluid dynamics [30], the vortex 
filament model is just a convenient idealization; the vorticity in a realistic 
classical fluid flow rarely takes the form of clearly discrete vorticity filaments. 
However, the vortex filament model is accurate and realistic for a quantized 
vortex in superfluid helium. 

3.2.1. Schwarz' s model 

In considering the velocity field created by a vortex filament in this sub- 
section, we introduce Schwarz's model, which is useful for superfluid helium. 

The incompressible velocity v(r) created by the vorticity source u>(r) 
satisfies divv = and rotv = u>, whose solution is [3o[ 

v « = ^ / ^ r ') x |7^ rfr ' ( 2 ) 

This representation is applied to the vortex filament formulation, which de- 
scribes a quantized vortex as a filament passing through the fluid, having a 
definite direction corresponding to its vorticity. The three-dimensional con- 
figuration of vortex filaments can be represented by differential geometry. A 
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point on a filament at a time t is represented by the parametric form s(q,t) 
with the one- dimensional coordinate q along the filament. The positive direc- 
tion of q is taken to be the direction of the vorticity. Then, ds/dc; = s'(q,t) 
is a unit tangential vector along the filament at s, and d 2 s/dq 2 = s"(q,t) is 
the principal normal vector at s with magnitude R~ x , where R is the local 
radius of curvature. Except for the thin core region, the superflow velocity 
field has a classically well-defined meaning and can be described by ideal 
fluid dynamics. The vorticity with quantized circulation k is focused only on 
the filament, represented by 

u(r,t) = K J s'{s,t)8{r-a{s,t))ik, (3) 

where the integration is taken along the filament. Inserting Eq. 03) into Eq. 
02]) yields the Biot-Savart expression 

k f s'k,t) x (r — sU,t)) , 
4vr Jc |r-s(c;,i)| 3 

Thus some configuration of vortex filaments gives the superfluid velocity field 
v s ,Bs(r,t). 

Considering the forces acting on the vortex filament, we derive the equa- 
tion of motion, namely Schwarz's equation. When a vortex filament moves in 
the superflow field v s , the effective forces are the Magnus force, the mutual 
friction force, and the inertial force. The Magnus force refers to the lift force 
acting on a spinning object when it moves in a fluid. The Magnus force for 
our vortex filament per unit length is written as 



p s Ks' x (s - v s ), (5) 



where s = ds/dt refers to the velocity of the filament at s. The Magnus force 
tends to move the vortex filament at s normal to both the vorticity ks' and 
the relative velocity s — v s between the vortex velocity and the superflow. At 
finite temperatures mutual friction works through the interaction between 
the normal flow and the vortex core: 



-ap s Ks' x [s ; x (v n - v s )] - a'p s Ks' x (v n - v s ), (6) 



where a and a' are coefficients dependent on temperature [28|. The right 
hand side represents the force normal to s'. We can write the equation of 
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motion of the vortex filament per unit length as 



d 2 s 

m &~Jp = *m + fo, (7) 

with the effective mass m e g- of the filament per unit length. The effective 
mass should be of the order of p s a,Q, which is usually quite small compared 
with other terms because of the small core radius ao- Thus we can neglect 
the inertia term, so that Eq. (J7J) is reduced to Schwarz's equation 

s = v s + as' x (v„ - v s ) - a's' x [s' x (v n - v s )]. (8) 

When we attempt to obtain v s at a point s(<To) along a filament from Eq. 
01]), the integral diverges as s(g) — > s(q ). In order to treat this difficulty we 



introduce the localized induction velocity proposed by Arms and Hama [31 
By using a cutoff R, the integral of Eq. 03]) is divided into the contribution 
within R around and that from the other distant region. The neighborhood 
of s(<j ) is represented by 

s(?) = B( ft ) + (? - %)s'( % ) + I(? - ft )V(<b). (9) 
This expression with Eq. (j3j) yields the local contribution 

Vi(ft) = £ r -B'(ft) X B ff (ft) = /3s'(, ) x s"( ?0 ) (10) 

with /3 = (k/47t) log(i?/ao). The parameter ao is the cutoff, corresponding 
to the core radius. This velocity is called the localized induction velocity or 
self-induced velocity. Since the contribution from the outer region is obtained 



by the usual integration, v s of Eq. (JH1) is reduced to [28 



n i ii K f ( s i — r) x dsi _ , . 

4vr Jc \si - if 

The second term represents the non-local field obtained by integrating the 
integral of Eq. (j3J) along the rest of the filament, except in the neighborhood 
of s. 

A better understanding of vortices in a real system is obtained when 
boundaries are included in the analysis. For this purpose, a boundary- 
induced velocity field v s b is added to v s , so that the superflow can satisfy the 
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boundary condition of an inviscid flow; that is, the normal component of the 
velocity should disappear at the boundaries. To allow for another, presently 
unspecified, applied field, we include v s>a . 

Consequently, the total velocity So of the vortex filament without dissi- 
pation is 

s = /3s' x s" + — / ^ + v S)6 s + v s , a s . 12 

4vr J c |si-r| 3 

At finite temperatures, it is necessary to take into account the mutual friction 
between the vortex core and the normal flow v n . Including this term, the 
velocity of s is given by 

s = s + as' x (v n - s ) - a's' x [s' x (v n - s )], (13) 

where So is calculated from Eq. ( Tl2l) . 

3.2.2. Basic motion of vortex filaments 

In this subsection we discuss the simple motion of vortex filaments in 
order to develop a basic understanding. The addressed equations of motion 
are Eqs. (TT2l) and (pi). 

The first term of Eq. ( fT2l) refers to the localized induction field arising 
from a curved line element acting on itself. The mutually perpendicular vec- 
tors s', s", and s' x s" are directed along the tangent, the principal normal, 
and the binormal, respectively, at the point s, and their respective magni- 
tudes are 1, Br 1 , and Rr 1 , where B is the local radius of curvature. Thus, 
the first term represents the tendency to move the local point s in the bi- 
normal direction with a velocity inversely proportional to B. Neglecting the 
non-local terms is referred to as the localized induction approximation (LIA). 
This approximation is believed to be effective for analyzing isotropic dense 
tangles due to cancellations between non-local contributions [29(. However, 
the LIA lacks the interaction between vortices, and is not necessarily suitable 
for the description of a realistic vortex tangle, as shown in 3.3. 

We consider vortex motion under mutual friction. The first example is 
the propagation of a vortex ring. Without any other vortices or a velocity 
field at zero temperature, a vortex ring with a radius B just propagates with 
a velocity approximately equal to the self-induced velocity (3/B normal to 
the circular plane, keeping its shape. We consider the motion of a vortex 
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ring under the applied fields v n and v Sj0 . Using the LIA and neglecting the 
friction term a' for simplicity, Eqs. ( TT2|) and (TT3"j) are reduced to 

s = (3s' x s" + v s , a + as' x (v n - v Sja - (3s' x s"). (14) 

The self-induced velocity is supposed to be in the z direction with s' x s" = 
z/R. If we take v„ = v n z and v s a = v s>a z, Eq. (JT4"j) describes the time 
development of R as 

a(v n — v aj a-^L\ (15) 



dt V" " u 

The absence of any applied fields gives dR/dt = —a((3/R), whose solution is 
R = \JRq — 2a(3t with the initial radius Rq. Thus the mutual friction shrinks 
the ring. The presence of fields complicates the situation. For v n — v s ^ a < 0, 
the ring always shrinks. When v n — v S)Cl > 0, the ring expands for v n — v s ^ a > 
(3/R and shrinks otherwise, which leads to a critical radius of curvature given 
by R c ~ (3 /{v n — v Sta ). These simple considerations emphasize the important 
role of the mutual friction. The mutual friction can both shrink and expand 
a ring depending on the applied field and the radius. In a vortex tangle, 
fine structure with a small radius of curvature generally shrinks under the 
mutual friction. In other words, the mutual friction tends to make the vortex 
configuration smooth. 

The second example we consider is the motion of two parallel or antiparal- 
lel vortices. Suppose that two straight vortices with circulation n are placed 
in parallel with a distance 2r. Each vortex moves through the velocity n/Airr 
from the other. At zero temperature, they rotate around their middle point, 
maintaining their distance. At a finite temperature, Eq. f fl3|) . neglecting the 
a' term, yields dr/dt = a(fi;/47rr), whose solution is r = yrjj + (aK/2n)t 
with the initial distance r . Thus the two vortices spiral outward, which 
means that their interaction is effectively repulsive. On the other hand, two 
antiparallel straight vortices move normal to the line segment connecting 
them with velocity n/Anr at zero temperature. The mutual friction at a 
finite temperature reduces their distance as r = a/^q — (aK/2ir)t. The in- 
teraction between them is effectively attractive, which eventually leads to 
pair-annihilation of two voritces. 

3.2.3. Numerical simulation 

The numerical simulation method based on this model has been described 



in detail elsewhere [28|, |29|, |32|, [33] . A vortex filament is represented by a sin- 



gle string of points separated by a distance A<^. The vortex configuration 
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at a given time determines the velocity field in the fluid, thus moving the 
vortex filaments according to Eqs. (IT2"|) and (1131) . When vortices move, they 
have chances to encounter other vortices. Thus, vortex reconnection should 
be properly included when simulating vortex dynamics. A numerical study 
of a classical fluid shows that the close interaction of two vortices leads to 



their reconnection, primarily because of viscous diffusion of the vorticity (34 
Schwarz assumed that two vortex filaments reconnect when they come within 
a critical distance of one another, and showed that statistical quantities such 
as the vortex line density were not sensitive to how these reconnections oc- 
cur 28|, |29| . Even after Schwarz's study, it remained unclear as to whether 
quantized vortices can actually reconnect. However, Koplik and Levine di- 
rectly solved the GP equation to show that two closely quantized vortices 
reconnect, even in an inviscid superfluid [35J. Therefore such an artificial 
procedure of vortex reconnection is currently thought to be allowed in the 
vortex filament model too. The more modern and reasonable procedure is 
to reconnect two vortices when they pass within the spatial resolution A^ 
with unit probability. Every vortex initially consists of a string of points at 
regular intervals of Ac;. When a point on a vortex approaches another point 
on another vortex more closely than the fixed space resolution A<j, we join 
these two points and reconnect the vortices. This reconnection procedure is 
standard in the vortex filament model, but a different procedure is used in 



some studies 36 



3.2.4- Some statistical quantities 

Some important quantities that are useful for characterizing the vortex 



tangle are introduced below 29|]. The vortex line density (VLD) is the total 



length of vortex lines per unit volume, defined by 



<k, (16) 



where the integral is performed over all vortices in the sample volume Q. The 
anisotropy of the vortex tangle that is formed under the counterflow v ns is 
represented by the dimensionless parameters 

h = ^7 f [l-(s'-rn) 2 ]*, (17) 
/^ = ^[l-(s'-r ± ) 2 ]*, (18) 
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ns 



Here, r\\ and represent unit vectors parallel and perpendicular to the v 
direction, respectively. Symmetry generally yields the relation I\\/2 + I± = 1. 
If the vortex tangle is isotropic, the averages of these parameters are 7ii = 
7j_ = 2/3 and 7/ = 0. At the other extreme, if the tangle consists entirely of 
curves lying in planes normal to v ns , then 7ii = 1 and I± = 1/2. 

When we perform numerical simulations, we should deal with such sta- 
tistical quantities as well as the dynamics of each vortex. The characteristic 
behavior of the VLD L in thermal counterflow was considered by Vinen. 
In order to describe amplification of a temperature difference at the ends 
of a capillary retaining thermal counterflow, Gorter and Mellink introduced 
some additional interactions (mutual friction) between the normal fluid and 
superfluid [16]. Through experimental studies of the second-sound attenua- 
tion, Vinen considered this Gorter-Mellink mutual friction in relation to the 



macroscopic dynamics of the vortex tangle [22|, [23|, |24], |25|. Assuming homo- 
geneous superfluid turbulence, Vinen obtained an equation for the evolution 
of L(t), which we call Vinen's equation: 

§ = ^|v ns |L^- X2 fL 2 , (20) 
at 2p 2tt 

where Xi is a constant, and B and X2 are temperature-dependent parame- 
ters. The first term represents the energy injection from the normal fluid 
to the vortices. The second term denotes the energy dissipation of vortices 
due to reconnection between vortices. The first and second terms indicate 
the growth and the degeneration of a vortex tangle, respectively. After the 
growth period of the VLD, the vortex tangle enters a statistically steady 
state. In the steady state, the VLD is obtained by setting dL/dt equal to 
zero, which gives 

L = iXs, (21) 
where 7 = ixBp n x\l KPX2 is a temperature-dependent parameter. This rela- 



tion can describe a large number of the observations of stationary cases [27 
When we conduct a simulation of the counterflow, the confirmation of Eq. 
f l2"Tj) is a touchstone. 

3.3. Thermal counterflow turbulence by the full Biot-Savart law 

The difficulty in accounting for the nonlinear and nonlocal dynamics of 
vortices has long delayed progress in achieving a microscopic understanding 
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of QT. It was Schwarz who made the breakthrough [29J. He investigated 
counterflow turbulence using the vortex filament model and dynamical scal- 
ing. The observable quantities obtained by his calculation agreed well with 
the experimental results for the steady state of vortex tangles. This study 
confirmed the idea proposed by Feynman that superfluid turbulence con- 
sists of a quantized vortex tangle. However, thermal counterflow turbulence 
was still far from being perfectly understood. The numerical simulation of 
Schwarz had serious defects. One is that the calculations were performed un- 
der the LIA neglecting interactions between vortices. Schwarz reported that 
as a result the layer structure is constructed gradually when periodic bound- 
ary conditions are applied. Of course, this behavior is not realistic. In order 
to address this, an unphysical, artificial mixing procedure was employed, in 
which half the vortices are randomly selected to be rotated by 90° around 
the axis defined by the flow velocity. This method enables the steady state 
to be sustained under periodic boundary conditions. These defects cause us 
to conjecture that the LIA is unsuitable due to the absence of interactions 
between vortices. 

Adachi et al. performed numerical simulations of counterflow turbulence 
using the full Biot-Savart law under periodic boundary conditions and suc- 
ceeded in obtaining a statistically steady state without any unphysical pro- 
cedures [iij]. Figure 1 shows a typical result of the time evolution of the 
vortices, whose VLD grows as shown in Fig. [2J The initial configuration 
consists of six vortex rings. 

In the first stage (0 < t < 0.4 s), the critical radius R c determines the vor- 
tex destiny. Vortex ring sections in which the radius of curvature exceeds R c 
expand in the direction perpendicular to v ns through mutual friction, while 
small vortex rings shrink. Thus, vortices evolve and become anisotropic. At 
the end of this stage, large vortices appear that are comparable to the system 
size under periodic boundary conditions. These vortices survive with a large 
radius of curvature, and continuously generate small vortices by reconnec- 
tions in the subsequent stages so that they function as "vortex mills" 37] . In 
the second stage (0.4 < t < 2.0 s), vortex tangles undergo continuous evolu- 
tion despite the decreasing anisotropy. As vortex rings expand, reconnections 
between vortices occur frequently. Reconnections generate vortices with var- 
ious curvatures, resulting in them shrinking and expanding as discussed for 
the first stage. Local sections with a small radius of curvature formed by re- 
connections have an almost isotropic self-induced velocity, which prevents the 
vortices from lying perpendicular to v ns . In addition, as the VLD increases, 
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Figure 1: Development of a vortex tangle by the full Biot-Savart calculation in a periodic 
box with a size of 0.1 cm. Here, the temperature is T — 1.9 K and the countcrflow velocity 
v ns — 0.572 cm/s is along the vertical axis, (a) t = s, (b) t = 0.05 s, (c) t — 0.5 s, (d) 
t = 1.0 s, (e) t — 3.0 s, (f) t — 4.0 s. [Adachi, Fujiyama and Tsubota: Phys. Rev. B 
81 (2010) 104511, reproduced with permission. Copyright 2010 the American Physical 
Society] 



vortex expansion becomes slower than in the first stage because the recon- 
nection distorts vortices, which prevents a vortex from smoothly expanding. 
In the third stage (t > 2.0 s), the statistically steady state is realized by the 
competition between the growth and decay of a vortex tangle. The growth 
mechanism is still vortex expansion through mutual friction. The decay 
mechanism either creates vortices with local radii of curvature smaller than 
R c or vortices with the self-induced velocity oriented in the opposite direction 
to v ns after the reconnections. The increasing VLD causes more reconnec- 
tions so that the decay mechanism becomes effective. When the VLD has 
increased sufficiently, the two mechanisms begin to compete so that the vor- 
tex tangle enters the statistically steady state. The LIA calculation cannot 
realize this competition, which shows that vortex interaction is essential for 
creating a steady state. 
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f(s) 

Figure 2: Vortex line density as a function of time for four different countcrflow velocities. 
[Adachi, Fujiyama and Tsubota: Phys. Rev. B 81 (2010) 104511, reproduced with 
permission. Copyright 2010 the American Physical Society] 



The obtained steady states almost satisfy the relation of Eq. (121]) when 
v ns and L are relatively large, as shown in Fig. [3j Table [1] shows the param- 
eter 7 as a function of T. The results quantitati vely ag ree with the typical 
experimental observations of Childers and Tough [271 . |38[ . Additionally, there 
is a critical velocity of turbulence, below which vortices disappear. This crit- 



it is 



ical velocity has been measured in many previous studies |27|, [39|, |40 
given by 

_ 2.5 + 1.44a 

Vns,c ~ ' 

where d is the channel size of the experimental system and a is a constant 
of order unity. In the simulation, the system size may be taken to be the 
size of the periodic box. Then, Eq. ( 122|) gives v nS)C ~ 0.1 cm, which is almost 
consistent with the numerical results. However, the temperature dependence 
of v nSjC should be considered. Equation (|22|) states that v ns>c should decrease 
with T, which differs from the behavior in Fig. EJ The numerical results show 
that v nStC decreases with T below 1.9 K but increases slightly at 2.1 K. This 
is because the strong mutual friction makes the vortices so anisotropic that 
they cannot form enough reconnections with other vortices, and so become 
degenerate. 

Thus a full account of the intervortex interaction by the full Biot-Savart 
law enables us to obtain the statistical steady states without any unphysical 
procedures. A detailed comparison between the numerical results of the full 



Biot-Savart law and the LIA is discussed in Ref. 33 
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Figure 3: Steady state vortex line density L(i) as a function of the counterflow velocity v ns . 
The error bars represent the standard deviation. [Adachi, Fujiyama and Tsubota: Phys. 
Rev. B 81 (2010) 104511, reproduced with permission. Copyright 2010 the American 
Physical Society] 



However, the situation may be not so simple. The above simulation 
was done under the assumption that the normal flow is laminar or uniform. 
By the recent visualization experiments using metastable helium molecules, 
Guo et al. showed that the normal fluid could be turbulent too at relatively 
large velocities (41] . In order to take account of the turbulent normal flow, 
we should address the coupled equations of the vortex dynamics and the 
Navier-Stokes equation describing the normal flow. It would be a future 
work. 



3.3.1. Velocity statistics 

Velocity statistics, namely the probability density function (PDF) of the 
velocity field, is another important statistic in turbulence. It is known that 
the PDF of classical viscous turbulence is Gaussian 
remains, what happens to the PDF in QT? 



42l . |43| . The question 
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T(K) 


7rmm(s/cm 2 ) 


7 exp (s/cm 2 ) 


'II 


1.3 


53.5 


59 


0.738 


1.6 


109.6 


93 


0.771 


1.9 


140.1 


133 


0.820 


2.1 


157.3 




0.901 



Table 1: Line density coefficients 7 and anisotropy parameter In . 7„ Mm and 
numerical results and the experimental results of Childers and Tough [27, S 
[Adachi, Fujiyama and Tsubota: Phys. Rev. B 81 (2010) 104511, reproduced with 
permission. Copyright 2010 the American Physical Society] 



, xp denote our 
, respectively. 



Paoletti et al. 44j performed visualization of quantized vortices in a re- 
laxation process of counterfiow using solid hydrogen particles and obtained 
the non-classical (non-Gaussian) PDF of the particle velocity. They reported 
that the non-classical statistics are due to the velocity induced by the recon- 
nection of a quantized vortex because the PDF exhibits a power-law dis- 
tribution of v ~ 3 , which is derived from the vortex velocity before or after 
reconnection. However, they observed the velocity of particles, which is not 
necessary the velocity of the superfiow. The non-classical velocity statistics 



were also confirmed by White et al |45j. They performed numerical simu- 



lations of QT in a trapped BEC by calculating the GP equation to obtain 
the PDFs of the superfiow field. The PDFs do not show classical Gaus- 
sian distributions, but rather power-law distributions, due to the velocity 
field v = K,/(2irr) induced by the singular quantized vortex, where r is the 
distance from the core of a quantized vortex. 

Adachi et al. studied the PDF of superfiow for the steady state of coun- 



terfiow |46J]. Figure H] shows the PDFs of v x , v y , and v z , where the counterfiow 
is applied along the z direction. The PDFs exhibit a non-Gaussian distri- 
bution with a large tail in the high-velocity region. Since the vortex tangle 
of steady counterfiow turbulence is isotropic in the direction perpendicular 
to the relative velocity v ns , the PDF of v x almost overlaps with that of v y , 
with the peaks of two PDFs at v x — and v y = 0. In contrast, since the 
superfluid velocity v sa = —0.496 cm/s due to counterfiow is applied in the 
—z direction, the PDF of v z has a peak at v sa . 

The PDF in the high-velocity region shows a power-law distribution of 
Pr(uj) oc v ~ 3 (i = x, y, z). For the single, straight vortex case, the probability 
of separation occurring between r and r + dr is 2irrdr, and the velocity 
scales as 1/r, which leads to Pr(t>) ~ Pr(r(v))\dr/dv \ ~ 1/v 3 . The PDF 
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Figure 4: (Color online) Probability distribution of the velocity components v x , v y , and 
v z in steady counterflow turbulence of v sa = —0.496 cm/s at T = 2. IK. The vertical 
dot-dashed line indicates v sa . [Adachi and Tsubota: Phys. Rev. B 83 (2011) 132503, 
reproduced with permission. Copyright 2011 the American Physical Society] 



converges to a Gaussian distribution in the low-velocity region, probably 
because the vortex configuration is random in the tangle. We can roughly 
estimate the transition velocity from the Gaussian distribution to the power- 
law distribution. In order to easily understand the velocity field induced 
by multiple vortices, we consider the simple case of two straight parallel 
vortices. Although the 1/r velocity primarily appears near each vortex, in the 
region halfway between vortices, the velocity becomes complicated because 
the velocities induced by the two vortices become comparable and interfere 
with each other. Hence, the statistics of velocity appear to change near 
the midpoint between vortices. In the vortex tangle, the mean inter-vortex 
distance is denoted by I = 1/L 1 / 2 , and so the midpoint between the vortices 
is located at 1/2. Thus, the transition velocity of the statistics should be 
represented by 

* = *m- (23) 

The transition in Fig. H] certainly occurs approximately at this scale. Thus 
the PDF shows the classical behavior at the low- velocity region and quantum 
behavior at the high-velocity region. 

3.4- Quantum turbulence created by vibrating structures 

Recently, vibrating structures, such as discs, spheres, grids, and wires, 



have been widely used for research into QT [47| . Despite detailed differences 
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between the structures considered, the experiments show some surprisingly 
common phenomena. 

This trend started with the pioneering observation of QT on an oscillating 



microsphere by Jager et al. [481] . The sphere used by Jager et al. had a 
radius of approximately 100 /im, and was made from a strongly ferromagnetic 
material with a very rough surface. The sphere was magnetically levitated 
in superfluid 4 He and its response with respect to the alternating drive was 
observed. At low drives, the velocity response v was proportional to the 
drive Fp, taking the "laminar" form Fn = X(T)v, with the temperature- 
dependent coefficient A(T). At high drives, the response changed to the 
"turbulent" form F D = j(T)(v 2 — Vq) above the critical velocity v . At 
relatively low temperatures the transition from laminar to turbulent response 
was accompanied by significant hysteresis. Subsequently, several groups have 
experimentally investigated the transition to turbulence in the superfluids 
4 He and 3 He-B by using^grids @ HuL H , wires @, E 
and tuning forks |60l 
a review article [4 




The details of the observations are described in 
Here we shall briefly describe a few important points 
necessary for the current article. 

These experimental studies reported some common behavior indepen- 
dent of the details of the structures, such as the type, shape, and surface 
roughness. The observed critical velocities are in the range from 1 mm/s 
to approximately 200 mm/s. Since the velocity is usually much lower than 
the Landau critical velocity of approximately 50 m/s, the transition to tur- 
bulence should come not from intrinsic nucleation of vortices but from the 
extension or amplification of remnant vortices. Such behavior is shown in 

Figure [5] shows 
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the numerical simulation by the vortex filament model 
how the remnant vortices that are initially attached to a sphere develop into 
turbulence under an oscillating flow. Such behavior must be related to the 
essence of the observations. 

Generally it is not easy to control the remnant vortices in an actual 
experimental setup. However, Goto et al. succeeded in preparing a vibrating 



wire free from remnant vortices [59]. This wire never causes a transition 
to turbulence by itself; it can cause turbulence only when it receives seed 
vortices from another wire. Such a simulation was performed by Fujiyama et 
al. as shown in Fig. [6] 63]. The sphere oscillates horizontally; the diameter 
of the sphere is 3 /im, the frequency of the oscillation is 1590 Hz, while the 
oscillation velocity is chosen in the range of 30-90 mm/s. Vortex rings of 
radius 1 /im are injected from the bottom of the medium [Fig. Efa)]. When 
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Figure 5: Evolution of the vortex line near a sphere of radius 100 /iin in an oscillating 
superflow of 150 mm -1 at 200 Hz. [Hanninen, Tsubota and Vinen: Phys. Rev. B 
75 (2007) 064502, reproduced with permission. Copyright 2007 the American Physical 
Society] 

the vortex rings collide with the sphere, reconnections occur and the vortices 
become attached to the sphere. Then, the attached vortices are stretched as 
the sphere moves [Fig. MJo) and (c)]. Due to the successive injection of vortex 
rings the process is repeated and the stretched vortices form a tangle around 
the sphere [Fig. IBfd)]. The vortices grow in size and some then detach from 
the sphere. In spite of the detachment, the oscillating sphere still sustains 
the vortex tangle when its velocity is relatively large. The vortex line length 
in a finite volume including the sphere was calculated at different oscillation 
velocities as a function of elapsed time. The loss of the vortices escaping from 
the volume balances the injection and the growth of the vortices so that the 
line length saturates. Only a slight increase in the line length can be observed 
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(c) t = 58 ms (d) t = 100 ms 



Figure 6: Time evolution of turbulence generation for the case of a sphere oscillating with 
a velocity magnitude of 90 mm/s. See the text for details. [Fujiyama and Tsubota: Phys. 
Rev. B 79 (2009) 094513, reproduced with permission. Copyright 2009, the American 
Physical Society] 
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for a velocity magnitude of 30 mm/s, which means that the vortices are not 
stretched by the sphere. The saturated value of the line length increases 
with the oscillation velocity magnitude. For velocity magnitudes above 50 
mm/s, the saturated line length value is much larger than the injection of 
vortices, which suggests that vortex tangles are formed around the sphere. 



This behavior is qualitatively consistent with the observations [59 



In order to characterize the transition to turbulence, Fujiyama et al. also 



studied the drag force 63]. The drag force acting on an object in a uniform 



flow is generally represented by 

I 

2 



d = ~C DP AU\ (24) 



where Cd is the drag coefficient, p is the fluid density, A is the projection 
area of the object normal to the flow, and U is the flow velocity. It is known 
in classical fluid mechanics how Cd depends on the properties of the flow. 
At low Reynolds number, Stokes's drag force acts on the object, which is 
proportional to the magnitude of U, with the result that is inversely pro- 
portional to U. When the flow becomes turbulent at high Reynolds number, 
Cd is of order unity. Fujiyama et al. estimated the drag force for the cases 
such as those in Fig. [6j The amplified line length can be related to the 
increase in energy, which should be equivalent to the work by the sphere. 
The drag coefficient Cd obtained by these considerations was of order unity. 
Thus we could confirm an analogy between CT and QT for this problem too. 

Another important simulation was performed for the current problem. 
Bradley et al. studied experimentally the transition to QT in the B phase 



of superfluid He [64J by vibrating a grid [52]. In superfluid He-B they set 



up a grid, 5.1 x 2.8 mm, composed of ~10 /mi square cross section wires 50 
pm apart. Directly in front of the grid were two vibrating wires, which could 
observe the vortices coming from the grid. The observed behavior showed 
two distinct regimes. At low grid velocities (below 3.5 mm/s) the two wires 
caught only vortex rings coming ballistically from the grid. At high grid 
velocities, however, the observation of two wires showed a signature of QT; 
the grid produced a vortex tan gle. Such behavior was confirmed through a 



simulation by Fujiyama et al. [65(. They followed the dynamics of vortex 
rings injected into a simulation "cell" such that the left-hand side of the cell 
represents the face of the grid. The simulation cell was a box of cross section 
200 /im x 200 fim and length 600 /im, as shown in Fig. [7J In the transverse 
directions the cell had periodic boundary conditions. Vortex rings of diameter 
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(a) t = 20 ms 



(b) t = 100 ms 




(c) t = 300 ms (d) t = 500 ms 

Figure 7: Simulation of QT formation. Each frame shows the vortex configuration at 
the labeled time. Rings injected from the left quickly collide and recombine to produce a 
vortex tangle. See the text for details. 



20 /im were injected at the left-hand side of the cell at a regular time interval 
Tj but at random positions and random angles within a ~20 deg. cone around 
the forward direction. The rings traveled at a self-induced velocity of 4.6 
mm/s. At a low injection rate (r» = 5ms) the simulation confirms that the 
rings travel essentially independently. At higher ring injection rates, however, 
corresponding to higher grid velocities, they found a very different behavior, 
as shown in Fig. [7] for n = 1.5 ms. Here the rings immediately start to 
collide and reconnect, establishing a vortex tangle, which corresponds to the 



behavior at high grid velocities in Ref. [52 



3.5. Gross-Pitaevskii model 

In a weakly interacting Bose system, the macroscopic wave function \1/ 
appears as the order parameter of Bose-Einstein condensation, obeying the 



Gross-Pitaevskii (GP) equation [66|, |67 



dt \ 2M y| 1 1 
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Here, g = 4irh 2 a/M represents the coupling constant characterized by the 
s-wave scattering length a, M is the particle mass, and p is the chemical 
potential. Expressing the order parameter with the amplitude and phase, i. e., 
= y/pe 1 ^ ', we obtain the condensate density p and the superfluid velocity 
v s = (h/M)V<f). The vorticity V x v s defined from the superfluid velocity 
vanishes everywhere in a singly-connected region of the order parameter, 
and all rotational flow is carried only by quantized vortices. The quantized 
vortex is defined as a topological excitation in which p vanishes at the core 
and rotates by 2ir around the core. The only characteristic length scale of 
the GP equation is the healing length defined by £ = h/ y/2Mgp with mean 
condensate density p; the vortex core size is given by £. 

The GP model can explain not only the dynamics of vortex lines but 
also vortex core phenomena such as reconnection and nucleation of vortices. 
However, strictly speaking, the GP equation is not applicable to superfluid 
helium, which is not a weakly interacting Bose system. The GP equation is 
well applicable to dilute atomic BECs 

3.5.1. Hydrodynamic properties and dynamics of a single quantum vortex 

Before investigating turbulent properties in the GP model, we briefly 
overview several hydrodynamic properties of the GP equation. With con- 
densate density p and superfluid velocity v s , the GP equation (1251) can be 
rewritten as follows: 

^ + V • (pv s ) = 0, (26) 



dt 2 Mp V 2 J 2M 2 V y/p )' 

Equations ( )26l) and (1271) express the equations of conservation of mass and 
momentum for a compressible inviscid fluid. The first term in the right- 
hand-side of Eq. f )27j) corresponds to an effective pressure p = gp 2 /2 due 
to the nonlinearity of the GP equation. The second term, the so called 
quantum pressure, has no analog in standard fluid mechanics, and becomes 
especially important at small scales comparable to the healing length such 
as near vortex cores, where p rapidly changes in the scale of £. Another big 
difference between QHD described by the GP model and perfect classical fluid 
hydrodynamics is the existence of quantized vortices. In the GP model, any 
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rotational flow is carried by quantized vortices with the quantized circulation: 



K = jv s -dl=— , (28) 

and it is well known that quantized vortices behave simply as vortex filaments 
in a perfect fluid. 



To study QT, we introduce a dissipation into the GP equation |68|, |69 
In atomic BECs, the main origin of the dissipation is considered to be in- 
teraction between the condensate cloud and the thermal component. In this 



section, we introduce the dissipation term in the following simple way [70 
We assume that the system is described by \l/ such that energy and particles 
are exchanged with a particle reservoir. The particle reservoir is a thermo- 
dynamic environment that lets the chemical potential of the system equal to 
that of the reservoir. The interaction with the particle reservoir is provided 
by an imaginary term in the GP equation: 

dt V 2M y| J K ' 

Here the wave function includes the chemical potential fi through the 
gauge transformation \& M = \I/e~* M *. Since T arises from the difference of the 
chemical potential, we can write T = 7(71 — /i pr ) with the chemical potential 
yU pr of the particle reservoir. We assume that the system is nearly in equilib- 
rium with the particle reservoir and that F is proportional to the difference 
in the chemical potentials. Using the approximation 

d^ 

~ (30) 

and n ~ /i pr , Eq. fl29l) becomes 

dV„ ( h 2 V 2 |2 



{t " 7) h m = ("W + 9 |vig + im ) ^' (31) 

Substituting ^ M = \l>e"* M ' into Eq. ( |3T|) . we finally obtain the modified GP 
equation with the effective dissipation 7: 

+ (32) 



26 



Introducing 7 conserves neither the energy nor the number of particles. For 
studying the hydrodynamics, however, it is sometimes realistic to assume 
that the number of particles is conserved. Hence, we can introduce the time 
dependence of the chemical potential so that the total number of particles 
N = J dr I ^ 1 2 is conserved. 

When 7<C 1, Eq. (j32J) becomes approximately 

if^r = (1 - n) + g |*| 2 - n) *, (33) 



dt v " V 2M 
which is widely known as complex-Ginzburg-Landau equation for the case 



when 7 is uniform [7l|. By using the density p and the superffuid velocity 



v s , Eq. (13"3"|) can be rewritten as 

| + V .(^)= 7 {»^-^-fcl!!^i, (34) 



^ 2 Mp\ 2 J 2M 2 \ J 2M { p 

The third term in Eq. ( )35l) works as a viscous term with the kinematic 
viscosity fry/(2M). We can therefore define the effective Reynolds number 

R Q t = (36) 

when we consider QT in the GP model. Here v s and 7 are the averaged 
values of v s and 7 over the space defined as 



dr \pv s \ / dr /ry 

, 7 = ^7 » (37) 



dr p dr p 



and Z) is the system size. 

Next, we investigate the dynamics of a single vortex line 68]. The steady 
solution of the GP equation with a single vortex line along the z axis is 
represented by 

*(r) = v^e^, (38) 



27 



where r = \J x 2 + y 2 and (p = tan 1 (y/x) are the radius and angle in cylin- 
drical coordinates. p{r) follows the equation 



h 2 



2M V dr 2 r dr 



VP 



9" 



hVp = o. 



(39) 



Starting from the solution of Eq. 
an external flow v P as 



h 2 



2M 



391) . we calculate the GP equation under 



(40) 



V 2 + g |*f -fit- ihv c ■ V * 



(41) 



and obtain the time development of the density p: 

d^fp d^fp 2 d^fp 

-WT = -Vcr—T- =F l V e<p + 7 V er ——, 

at dr dr 

which is second order in 7. Here v e = v er r + v etp cp. The vortex moves from 
the positive to negative sides of dy/p/dt in Eq. (14T|) . Assuming Jp~ oc r 
around the core, we obtain the time development of the vortex position r§. 



dr 

—— = v e ± 7i> c x 2 
dt 



(42) 



The first term states that the vortex moves with the velocity field v e |72| . The 



second and third terms describe the drag forces perpendicular and parallel 



to v r 



73, 74 



Recognizing v e in Eq. ( 142]) as the velocity field induced by other vortices, 
we can derive the effective vortex dynamics for a system with many vortices. 
In that case, the second and third terms work as the mutual friction forces 
by comparing Eq 



20, 21, 28 



with Eq. OH]) of the vortex-filament model 
with the corresponding coefficients a = 7 and a' = 7 2 . 

Furthermore, there is another important vortex dynamics which is not di 



35, 75, 76, 77 



rectly described in Eq. (|4*2]) : the reconnection of two vortices 
When two vortices are close to each other, they approach and become locally 
anti-parallel, and then reconnection occurs. Figure [S] shows an example of 
vortex reconnection given by the numerical calculation of Eq. (j25p starting 
from two straight vortex lines in a skewed orientation. Reconnection occurs 
even for the original GP equation (125]) without dissipation. We show that 
this process does not violate Kelvin's circulation theorem; moreover, the pro- 
cess becomes irreversible due to the emission of compressible excitations that 
have a wavelength smaller than £ 77 
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(a) (b) (c) (d) 




Figure 8: Reconnection of two vortices starting from a skewed position of two straight 
vortex lines. Simulation of the GP model of Eq. ([251) . (a) Initial state, (b) State just before 
the two vortices connect, (c) Connection of two vortex lines, (d) State after the vortex 
lines separate in a newly connected arrangement, also called reconnection. The contours 
in all the figures show the point of low density (10% of maximum density). [Kobayashi 
and Tsubota: J. Phys. Soc. Jpn. 74 (2005) 3248, reproduced with permission. Copyright 
2005 the Physical Society of Japan.] 



3.5.2. Quantum turbulence at zero temperature 

In this section, we overview our study of QT at zero temperature through 
the analysis of the GP equation. QT is defined as the turbulent state of a 
quantum fluid with highly tangled quantized vortices, and the term is often 
used to emphasize that the state is dominated by the behavior of quantized 
vortices at very low temperatures at which the thermal effect is negligible 

Q. 

For the analysis of turbulence, the energy spectrum is one of the most 
important statistical quantities. The energy spectrum can be obtained from 
the Fourier transformation of the equal-time two-point velocity correlation 
function: 

F(k, t) = ^Jdr e ik r J dr' (v{r', t) ■ v(r + r', t)) . (43) 

Here k is the wavenumber from the Fourier transformation and (• • ■ ) shows 
the ensemble average over statistically equivalent states. When studying 
the energy spectrum of QT, v is regarded as the superfluid velocity v s = 
(h/M)V(f). The integration of F(k,t) over the angle in wavenumber space is 
defined as the energy spectrum: 

E{k,t) = j^y J d<p k d6 k k 2 sm6 k F{k,t), (44) 

where 8 k and refer to the polar and the azimuthal angles in wave number 
space. The energy spectrum holds the following relation with the spatial 
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integration of the kinetic energy: 

J dkE(k,t) = -^ j 'dk (\v\ 2 ) = \ J dr (\v\ 2 )=E(t). (45) 

Here E(t) is the total kinetic energy per unit mass and v is the Fourier 
transformation of v. 

For turbulence of a typical viscous fluid, referred to as classical turbulence 
(CT) in this section, Kolmogorov proposed a statistically steady state of fully 
developed turbulence [78|, |79| in which energy is injected into the fluid at 
scales comparable to the system size D in the energy-containing range. In 
the inertial range, this energy is transferred to smaller scales without being 
dissipated, supporting a statistical law for the energy spectrum known as the 
Kolmogorov law: 

E(k) = Ce 2/3 k~ 5/3 . (46) 

The energy transferred to smaller scales in the energy-dissipated range is 
eventually dissipated at the Kolmogorov wavenumber kx = (e/u) 1 ^ through 
the kinematic viscosity v of the fluid, at the dissipation rate e. The Kol- 
mogorov constant C is a dimensionless parameter of order unity. 

The inertial range is thought to be sustained by the self-similar Richard- 
son cascade in which large eddies are broken up into smaller ones [80]. In 
CT, however, the Richardson cascade is not completely understood because 
it is impossible to definitively identify each individual eddy. In contrast, 
quantized vortices in QT are definite and stable topological defects. There- 
fore, QT gives the real Richardson cascade of definite quantized vortices, and 
thus is an ideal prototype for studying statistics such as the Kolmogorov law 
and the Richardson cascade in the inertial range of turbulence. Quantized 
vortices at finite temperatures can decay through mutual friction with the 
normal fluid at any scale. At very low temperatures, on the other hand, vor- 
tices can decay by the emission of compressible excitations and Kelvin waves 
through vortex reconnections. Therefore, dissipation occurs only at small 
scales; for large scales, we can obtain the turbulent state at high Reynolds 
number. 

We now consider the energy spectrum of the GP equation. The total 
energy of the GP equation per unit mass is 

E = — [dr ( — | W| 2 + . (47) 



MN / \ 2M 1 2' 
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Here N = J dr 2 is the total number of particles. E can be separated into 
the gradient energy -E gra d and the interaction energy E int . E gTlld is further 
separated into the kinetic energy E kin and the quantum energy E q as in the 
following: 



= E M „ +Et = JL.Jdr |V«f , E mt = ^ fir |* 

i r h 2 f 

E kin = dr p 2 , E Q = — / dr (V,fp) 2 , 

2N J q 2M 2 N J v VP; 



4 

(48) 



with p = ^v s . Efo n can be further divided into a compressible part E^ in due 
to compressible excitations and an incompressible part E^ in due to quantized 



vortices 81, 82 



E t=^Jdr {[ P r) 2 . (49) 

Here [•■■]° denotes the compressible part, i.e., V x [•••] c = 0, and 
denotes the incompressible part, i.e., V ■ [■ • - ] c = 0. The compressible part 
A c and the incompressible part A 1 of an arbitrary vector field A are given by 

A C = J2 ^T keikr i A [ = A- A c , (50) 

k 

where A is the Fourier component of A. Corresponding to each energy, 
there are several kinds of energy spectra. The most important is the energy 
spectrum of the incompressible kinetic energy: 

EUk ' t) = 2(2^N J d ^ d ^k 2 sm9 k (pf, (51) 

because it should obey the Kolmogorov law with the Richardson cascade of 
quantized vortices. Here p 1 is the Fourier transformation of the incompress- 
ible momentum density [p]\ 

Starting from a Taylor-Green vortex, Nore et al. simulated decaying 



turbulence by numerically solving the GP equation (1251) }8lL |82| . After some 
time, the initial vortices became tangled and the calculated energy spectrum 
El in obeyed the power- law behavior E^ n (k,t) oc k~' n . When vortices formed 
a tangle, the exponent r\ was about 5/3, but this value did not hold for long 
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because the turbulence was decaying with the conservation of total energy E 
in Eq. (H71) . Because of energy conservation, the energy of the vortices E l khi 
was transferred to compressible excitations E^ in with wavelength comparable 
to £ through repeated reconnections. Therefore, the dynamics of QT are 
affected by many compressible excitations and we cannot see the proper 
dynamics of quantized vortices, which causes El in (k,t) to deviate from the 
Kolmogorov law. 

To obtain QT free from compressible excitations, we use the modified 



GP equation ( 132]) discussed in the previous section [68|, |69[. The space- 
independent constant 7, however, acts on vortices as the mutual friction at 
finite temperatures, as discussed in Eq. (112"]) . and is inappropriate to study 
QT at zero temperature. Thus, we now consider the space dependence of 7 
as follows. The Fourier transformation of Eq. (132]) is 



i-7)fi— = [——-n )ty + gY. (52) 



dt \2M 

Here 7 and Y are defined to satisfy the relations: 

dr e ~ ife r 7^ = 7^-, I dr e- ik r m 2 m = Y, (53) 
dt dt J 11 v ; 

where both depend on Especially, 7(fe) directly acts on and dissipates the 
fe-component of \& '. When we choose 7 as the step-function form: 

7 = 7o0(*-27r/O, (54) 

we can expect that only compressible short-wavelength excitations produced 
via reconnections are dissipated. Because the system is dissipationless at 
scales exceeding £, we can investigate proper vortex dynamics at zero tem- 
perature in large scales without dissipation. We investigated the vortex dy- 
namics in the modified GP equation (152]) with the step function type dissipa- 
tion ([51]) and found that this kind of dissipation does not work as the mutual 
friction as described in Eq. ( I42p and just removes the short- wavelength com- 
pressible excitations. 

We next discuss simulations of QT described by Eq. f ]52]) with the dis- 
sipation of Eq. (]53]1 . Here we consider two kinds of turbulence: decaying 



turbulence 68] and steady turbulence [69]. For decaying turbulence, the ini- 



tial configuration was set to have a uniform density p = 1 and the phase 
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had a random spatial distribution. The random phase is generated by 
placing random numbers between — tt to 7r at every distance A and connect- 
ing them smoothly, which represents an energy injection at the scale of A. 
Because the initial superfluid velocity v s = (h/m) V0 given by the initial ran- 
dom phase is random, the initial wave function is dynamically unstable and 
soon produces turbulence with many vortices (see Fig. [5]). We confirm that 
only the compressible kinetic energy E^ in is decreased by the dissipation term 
and that the incompressible kinetic energy E^ in dominates the total kinetic 
energy E\^ n , demonstrating that only compressible excitations are effectively 
dissipated by the dissipation term 7. In the middle stage of the decay, the 
dissipation rate of the incompressible kinetic energy e^ in = —dE^ n /dt takes 
an almost constant value, showing the quasi-steady state of QT. In this range, 
the energy spectrum El in (k,t) is consistent with the Kolmogorov law: 

^in(^,t) = C(4 m ) 2/3 ^ 5/3 , (55) 

with the Kolmogorov constant C = 0.32, and smaller than that in CT which 
is estimated to be 1.4 < C < 1.8. Araki et al. obtained a Kolmogorov 
constant C ~ 0.7 in their numerical simulation of the vortex- filament model 
and this is also smaller than that in CT. This small Kolmogorov constant 
may therefore be characteristic of QT [83]. 

We also considered steady QT, which can be modeled by introducing en- 
ergy injection to the system. The advantages of steady QT over decaying 
QT are the following. First, steady turbulence gives a clearer correspondence 
with the Kolmogorov law: this is because the original statistics have been 
developed for steady turbulence. Second, it enables us to confirm the pres- 
ence of the energy- containing range, the inertial and the energy-dissipative 
range of QT. Third, in all ranges we can obtain the time-independent energy 
flux in wavenumber space. Consequently, it is possible to reveal the cascade 
process of QT, which occurs by quantized vortices. 

For energy injection at large scales, we introduce the moving random 
potential V(r,t) in the GP equation: 

(i -^" l 8T = (-W + sl * 1 2 -" + ^)*- < 56 > 

or its Fourier-transformed form: 

( i -^9T=(w-")* + 9i ' + K (57) 
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Figure 9: Snapshot of a vortex configuration in the turbulent state. The simulation box 
was set to (32£) 3 . Visualization of quantized vortices can be done using the following 
method. For four numerical grid points (x y z), (x + Ax y z), (x + Ax y + Ay z), and 
(x y + Ay z), we can calculate the phase shift. This value takes the values 0, 2tt, or — 2tt, 
with the two latter values indicating that a vortex line pierces this plaquette. Therefore, 
by calculating the phase shift at six plaquettes for each unit cube and taking the isosurface 
of its absolute value, we can visualize vortices. [Kobayashi and Tsubota: J. Phys. Soc. 
Jpn. 74 (2005) 3248, reproduced with permission. Copyright 2005 the Physical Society 
of Japan.] 



Here V is defined to satisfy the relation: 

dr e~ ik r V^ = V. 



(58) 



We set the statistical properties of V(r,t) to obey the Gaussian two-point 
correlation: 



(V(r,t)V(r',t'))=V 2 exp 



[x 



x 



2X1 



it - if) 
2T 2 



(59) 



This moving random potential has the characteristic spatial scale X and thus 
quantized vortices of scale X are nucleated when V is strong and T is short 
enough. We define the wavenumber separating the energy- containing range 
and the inertial range as 2k/Xq. The wavenumber 27r/£ between the inertial 
range and the energy-dissipative range is defined by the dissipation term 7. 
Therefore, our steady QT has an energy-containing range of k < 2ir/Xo, in- 
ertial range of 2tt/X < k < 27r/£, and energy-dissipative range of 27r/£ < k 
(see Fig. [TU]) . Starting from the uniform state p — 1 and <fi = 0, we develop 
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Figure 10: Flow of the incompressible kinetic energy EL (upper half of diagram) and com- 
pressible kinetic energy (lower half) in wavenumber space. [Kobayashi and Tsubota: 
J. Phys. Soc. Jpn. 74 (2005) 3248, reproduced with permission. Copyright 2005 the 
Physical Society of Japan.] 



the GP equation f )57|) with potential (159]) . and obtain steady QT in which 
E, E q , E^ in , and Ey dn are statistically steady. By choosing the appropriate 
parameters, the incompressible kinetic energy El in is always dominant in the 
total kinetic energy E^i n ; the introduced potential contributes to the nucle- 
ation of vortices rather than that of compressible excitation with long wave- 
length. The obtained energy spectrum was consistent with the Kolmogorov 
law in the inertial range. 

We further calculate two other important values: the energy dissipation 
rate £^ n and the flux of the incompressible kinetic energy from small 
to large wavenumbers. £^ n in steady turbulence can be equated to £y n = 
—dE^/dt after switching off the moving random potential. This is because 
the incompressible kinetic energy E^ in decays to the energy of compressible 
short-wavelength excitations. On the other hand, the energy flux II can be 
calculated by considering the scale-by- scale energy budget equation, which 
can be obtained by the time development of the cumulative incompressible 
kinetic energy: 

4n=2^ / dr CW) 2 - (60) 
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Here L k is the operator for the low-pass filter: 



LkW)] 



(27T) 



k<k 



dk / dr' e ih {r ~ r,) s(r) 



The time derivative of £{. in gives the energy budget equation: 



kin 



<9f 



+ n: 



kin 



kin + ^kin ^\in- 



kin 



kin' 



Here we introduce the cumulative energy injection J^ in 

1 



T 1 

J kin 



MN 

the cumulative energy transfer 7^ in 



dr L k \p] x ■ L k [^Vl/f 



7Sn=^ / drLfclpp-L, 



M \2M y /p 



99 



v s v s ■ Vp 



the cumulative energy dissipation PL 



kin 



Pi 



kin 



MiV 



+ 7^ 



2 



V 2 ^ M 2 v 



7 



/i 2 



£ 2 



+ 0( 7 2 



and the energy flux U\ 



kin 



nU = — I drL k [pf. L k [pV • ^ s + VpV«s] 1 . 



(61) 



(62) 



(63) 



(64) 



(65) 



(66) 



Equation (I62j) can be interpreted as follows: at a given scale k, the rate of 
change of the incompressible kinetic energy is equal to the energy injection 
by the force J^ in plus the energy transfer between vortices and density fluc- 
tuations 7^ in minus the energy dissipation P^ in minus the energy flux nj cin to 
smaller scales. 

Figures [U] (a) and (b) shows the energy flux nj^ and the energy dissipa- 
tion rate e\, in , and the energy spectrum E^ n (k) respectively for numerically 



36 



obtained steady QT. The energy flux lTj^ is nearly constant and is consis- 
tent with the energy dissipation rate e l kin in the inertial range, which indicates 
that the incompressible kinetic energy steadily flows in wavenumber space 
through the Richardson cascade at the constant rate H kin , and finally dissi- 
pates to compressible excitations at the rate £y n ~ n^ in . This energy flow is 
shown in the diagram of Fig. [10J The energy spectrum E kin (k) shown in Fig. 
[TT1 (b) is quantitatively consistent with the Kolmogorov law in the inertial 
range. The resulting Kolmogorov constant is C ~ 0.55, which is smaller than 
that in CT as well as for decaying turbulence. 
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Figure 11: (a) Dependence of the energy flux IlJ !in on the wavenumber k and the energy 
dissipation rate obtained from e kin = —dE kin /dt. (b) Energy spectrum E kin (k). The solid 
line is the Kolmogorov law C(e kin ) 2 / 3 k~ 5 / 3 . A simulation was performed in a periodic 
box of size 32 with parameters 70 = 1, Vq = 50, X — 4, and T — 6.4 x 10~ 2 . Here, 
length, wavenumber, energy, and time are normalized by £, l/£, h 2 / (2M£ 2 ), and (2M^ 2 )/fi 
respectively. With these parameters, the system enters steady turbulence after t ~ 25. 
nicin' e Lin ano - ^kin(^) were obtained from an ensemble average of 50 randomly selected 
states at t > 25. [Kobayashi and Tsubota: J. Phys. Soc. Jpn. 74 (2005) 3248, reproduced 
with permission. Copyright 2005 the Physical Society of Japan.] 



3.5.3. Quantum turbulence at finite temperatures 

Experimental studies of superfluid 4 He measured the energy spectrum 
of QT at finite temperatures, and suppor ted the Kolmogorov spectrum di- 



rectly or indirectly [84j, |85|, |86|, |87|, |88|, |89|, |90| . These experiments were also 



consistent in that they showed similarities between QT and CT. Vinen theo- 
retically considered this similarity and proposed that the superfluid and the 
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normal fluid are likely to be coupled by mutual friction at scales larger than 



the intervortex spacing I and would thus behave like a classical fluid [91 
This idea was confirmed by Kivotides et al. through numerical simulation of 
coupled dynamics of a vortex filament and a normal fluid |9~2| and by L'vov 
et al. through theoretical analysis of the two-fluid model [93J]. 

Using the GP model, the easiest way to discuss QT at finite temperature 
is to solve the modified GP equation (j3"2l with a constant 7, because, as we 
discussed earlier, the constant 7 behaves like the mutual friction with relation 
to the coefficients 7 = and 7 2 = a' (see Eq. (142]) ). As other methods to 



simulate atomic BECs at finite temperatures, the projected GP equation [94, 



95| . the stochastic GP equation [96(, and a coupled formalism involving the 
GP equation for condensate atoms and the Boltzmann equation for thermal 
cloud [97| are proposed. All methods include not only dissipation but also 
thermal fluctuation, which becomes more important at high temperatures. 

Before investigating QT at finite temperatures by using the above meth- 
ods, we have to clarify the origin of the dissipation term 7 in Eq. (132|) from 
the microscopic point of view and its temperature dependence [98|. In quan- 
tum fluids, dissipation comes from the interaction between the condensate 
and its excitation such as quantum and thermal fluctuations. For atomic 
BECs, the dynamics of the condensate and excitations can be describ ed by 



the GP and Bogoliubov-de Gennes (BdG) equations, respectively, [99|, llOO 
Our goal in this section is to microscopically clarify the dissipation mecha- 
nism of a quantum fluid with quantized vortices by numerically solving the 
coupled equations involving the GP and BdG equations. To do this, we start 
from the many-body Hamiltonian for bosons: 

H = fdr& [-~ - 11 + *. (67) 

Here ^ is the boson field operator. The time development of \& can be 
described by 

^fi 9^ = (-w-'' +9 * ,,1 ')*• (68) 

In the Bose condensed sy stem , the field operator \l> can be separated in terms 



of the mean-field ansatz lOll . 1 1021 

* = + x + C- (69) 
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In addition to the macroscopic wave function * = O{y/N /V), we define 
the first-order excitations x = 0(1/ yV) and the higher-order excitations 
( = 0(1/ y/N V) with the number of condensate particles N and the volume 
of the system V. Substituting Eq. (1B9l into Eq. (1681) and neglecting (, we 
obtain the GP equation: 



dt 



h 2 V 2 



2M 



(70) 



and the BdG equation: 

lh m = 



n 2 v 2 

' 2M 



(71) 



When the GP equation (FTP]) is expressed as ihd^/dt = Hqp^/, the corre- 
sponding Hamiltonian Hqp has the following imaginary term: 



lm[H t 



GP 



= —7 = Im 



(X 2 ) V 



(72) 



This defines the dissipation 7 of the condensate caused by the interaction with 
the noncondensed particles. We can calculate the dissipation 7 from Eq. ( 1721 
by numerically solving the coupled GP (170|) and BdG (17T]) equations. The 
BdG equation (1711 can be solved by using the Bogoliubov transformation: 



X 



Yl \ U J & 3 + r 'i'\, 



(73) 



where Uj and Vj are the Bogoliubov coefficients and dj and at are the an- 
nihilation and creation operators of a quasiparticle, respectively, for the jth 
energy level. Here, we assume that the quasiparticles are coupled with a heat 
bath at temperature T and they are in a locally equilibrium state: 



djdj 



e EilT _ X 



N 



(74) 



with the excitation spectrum Ej of quasiparticles. From these assumptions, 
we can expect that the energy of vortices or compressible excitations formed 
in the condensate is transferred to quasiparticles and finally dissipated to the 
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heat bath. Using Eqs. (l73p and f[T4l . we can obtain the final form of the 
coupled GP and BdG equations: 



-u + g(\^\ 2 + 2n e ) } + gm e m* 





h 2 v 2 

2M 


ih 9Uj ~ 
dt 


( h 2 v 2 

V 2M 


dv L = _ 


f h 2 v 2 


dt 


{ 2M 



- fi+ 2g\^\ 2 uj - g^ 2 Vj = A 



j- 



fi + 2g\^\ 2 jv j + g(^*) 2 u j = B j , 
{\uj?Nj + \vj\ 2 {Nj + 1)} , m c = - {ujV*{2Nj + 1)} 



(75) 



Ej = J dr Re (u*Aj + v*B 3 ) , 7 = Im 
For a given initial condition of we adopt the uniform excitations given by 

J(k r r) : 



Uj = e 



1 h 2 k 2 /(2M) + g\ty\ 



2V 



+ 1, 



Vj = e 



-i(fcj-r) 



1 h 2 k 2 /{2M)+g\m\ 



(76) 



2V 



E; 



1, 



as the initial condition for Uj and Vj. Here kj = 2ir(j x ,j y ,j z )/\/V, with 
nonzero integers j x , j y and j z . 

First, we attempt to calculate the coefficients of mutual friction as func- 
tions of temperature 0, 21, 28 1. When one straight vortex along the z axis 
is placed under the velocity field v e = (v e ,0, 0), the dynamics of the vortex 
position s(t) = (s x (t),s y (t),0) are described through Eq. (jS]) by 



s(t) = (s x (0) + (1 - a')v e , s y (0) - av e , 0) 



(77) 



Starting from the state \I/ in Eq. (I38p with one straight vortex line, we 
numerically solve the coupled GP and BdG equations under the velocity 
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field v c : 



cM> f h 2 V 2 1 
ih— = I —^f - A* + ffd^l 2 + 2n c) ~ iv e ■ V ^ + gm e m\ 

du- ( h 2 V 2 \ 
ih-^- = (- — - y. + 2g\^\ 2 - iv e ■ VJ u 3 - g* 2 v jt (78) 

dv- ( h 2 W 2 \ 
ih-^- = - I --^j- -11 + 2g\m\ 2 - iv e ■ V vj + g(V*) 2 Uj. 

We can obtain a and a' by comparing the position of the vortex with Eq. 
(177)) and find their monotonic increase with temperature as shown in Fig. 
H2J which is qualitatively consistent with mutual friction in superfluid 4 He 
at temperatures much lower than the superfluid critical temperature. This 
temperature dependence of a and a' may be a standard scale for measuring 
the temperature in atomic BECs with quantized vortices. 




Figure 12: Temperature dependence of the mutual friction coefficients a and a' . The plots 
represent the numerical results and the lines indicate fitting. A simulation was performed 
in a periodic box of size 4 with the velocity field v c — 0.1, where length, velocity, and 
temperature are normalized by £, H/(M£), and the critical temperature for the BEC for 
free bosons, respectively. [Kobayashi and Tsubota: Phys. Rev. Lett. 97 (2006) 145301, 
reproduced with permission. Copyright 2006 the American Physical Society] 

Next, we calculate the dissipation term 7 for the turbulent state. We 
begin with \& that includes several randomly placed vortices, as shown in 
Fig. [13] (a). After some time evolution, we calculate the dissipation term 
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7 in Eq. (T75|) . Figure d3] (b) shows the Fourier-transformed dissipation 
7 at several temperatures. At low temperature, dissipation works only at 
wavenumbers greater than 27r/£, which is consistent with the dissipation term 
j(k) = Jo0(k — 2n /£) used in the previous section for the simulation of QT 
at zero temperature. From this result, we expect that only short-wavelength 
excitations emitted during vortex reconnections or by high frequency Kelvin 
waves become dissipated at scales smaller than £. On the other hand, as 
the temperature increases, dissipation works at small wavenumbers as well, 
which is consistent with the above simulation of a single vortex, because 
dissipation at small wavenumbers acts as the mutual friction, as discussed in 
Eq. 



(a) 




0.2 



0.15 



& 0.1 J 



(b) 



0.05 





— T 
--T 

T 


= 0.01 
= 0.1 
= 0.5 




















5 10 


15 


20 



25 



Figure 13: (a) Example of the configurations of quantized vortices at t = 0. (b) Wave 
number dependence of the Fourier transformed dissipation term 7 at t — 1. A simulation 
was performed in a periodic box of size 4. Length, wavenumber, time, and tempera- 
ture are normalized by £, (2M£ 2 )/ft, and the critical temperature for the BEC for 
free bosons, respectively. 7 is obtained from an ensemble average of 25 different initial 
states. [Kobayashi and Tsubota: Phys. Rev. Lett. 97 (2006) 145301, reproduced with 
permission. Copyright 2006 the American Physical Society.] 



3.5.4- Two-dimensional turbulence 

In the previous section, we numerically verified that QT shows similar or 
the same statistical properties characterized by the Kolmogorov law. The 
physical explanation for this similarity is that quantization of vortices is not 
essential at scales larger than the mean intervortex spacing /, and vortices 
behave like eddies in CT forming vortex bundle structures. Considering this 
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situation, a new question arises: what happens during QT in a 2D system, 
in which quantized vortices take a point structure. 

Before discussing 2D QT, we consider 2D turbulence of a classical incom- 
pressible fluid obeying the Navier-Stokes equations 

v-« = o, 

9v , 1„ „ 2 (79) 

where p is the pressure and v is the kinematic viscosity. In the dissipationless 
limit v — > 0, the kinetic energy 

E= l -Jdrv\ (80) 

becomes motion invariant, and this plays a key role in the energy cascade 
because viscosity does not apply and the energy is not dissipated but is 
transferred from larger to smaller scales. In 2D system, the enstrophy Q is 
defined as 

Q= [ dr (V x vf (81) 



and also becomes motion invariant in the dissipationless limit, giving fur- 
ther cascade physics. Kraichnan recognized that the mot ion invariance of 



fl drastically modifies the physics of 2D turbulence [103] . There are two 
inertial ranges, the first for the cascade of the kinetic energy and the sec- 
ond for the enstrophy. In the energy cascading region, the direction of the 
energy flux is different from that in 3D turbulence, i.e., an energy cascade 
from small to large scales. The energy spectrum takes E(k) oc e 

2/3 fc -5/3 

and E{k) oc i] 3 ^ 2 k~ 3 (plus certain logarithmic corrections, unimportant in 
the current discussion) in the energy and enstrophy cascading regions, re- 
spectively. If energy is injected into a fluid in a wavenumber scale of ki, the 
inertial ranges for the energy cascade and the enstrophy cascade are formed 
in the wavenumber regions of k < ki and ki < k < k u , respectively, as the 
steady turbulent state, where k v is the viscosity cutoff. These predictions 
have been confirmed in laboratory experiments and using large-scale direct 



numerical simulations of Eq. f J79|) 104j . 1 1051 . Il06 



Our question is whether the inverse energy cascade and the enstrophy 
cascade are features of 2D QT described in the 2D GP equation. The qual- 
itative answer is that this is not necessarily the case. In a Bose condensed 
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system described by the GP equation, the enstrophy Q coincides (up to a 
prefactor) with the total number of vortex points. In the GP dynamics, how- 
ever, the total number of vortices is not conserved because vortices appear 
and disappear through pair creation and annihilation. Therefore, Kraich- 
nan's arguments become, generally speaking, irrelevant for 2D QT. Thus, for 
a certain range of parameters, where the creation or annihilation of vortex 
pairs becomes dynamically important, we can expect a direct energy cascade 
in 2D QT, exactly as in 3D turbulence. The same argument can be made 
based on the Euler equations. As shown in Eq. (|27|) . the GP equation is 
reduced to a compressible Euler equation with an effective pressure func- 
tion P = gp — /? 2 (V 2 v /p)/(2M v /p). However, a compressible Euler equation 
does not preserve the enstrophy. Therefore, at some level of compressibility 
(characterized by the Mach number, the ratio of the turbulent velocity fluc- 
tuations to the sound velocity), the direction of the energy flux can change 
its sign and, instead of an inverse energy cascade, we expect a direct cascade 
typical for 3D turbulence. 

Experimentally, we can assume a 2D BEC that is dynamically frozen in 
the transversal direction, i.e., a pancake-shaped BEC that is strongly trapped 
along the z axis. Let us consider a trapped BEC in the harmonic oscillator 
potential described by the GP equation 




Here, w x ,y,z are the oscillator frequencies on the x, y, and z axes. For sim- 
plicity, we here assume u x = u y = uj_. The strength of the trap along the z 
axis is determined by the ratio u z /u±, and a pancake-shaped potential gives 
oj z /u>± ^> 1. If f\jjj z is much larger than the chemical potential /i, the wave 
function can be separated into *&±(x, y) in the x-y plane and ^ z {z) along 
the z axis approximated by a Gaussian function as 

e -*7(2«2) 

m = tf z tf ± = qr ± . (83) 

Here, a z = \J hj (Moo z ) is the trap length along the z direction. Substituting 
Eq. (E3) into Eq. (E2D, we obtain 

. f h 2 W\ M. 2 2 2 2 , |T , 2 1 T ,o,N 

= [~^f ~ »± + Y^- x + ^ + 9±\*±\ 2 > *±- (84) 
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Here, = d 2 /dx 2 + d 2 /dy 2 , p± = p — ftw z , and g± = gj{\p2-na z ). As in 
the 3D case, we can define the 2D density p± = \^±\ 2 and 2D healing length 
a = fi/V2Mg ± p ± . 



<»-sr = -^nr - *l + »±l*±l a *±- (85) 



For simplicity, in the work |107j Numasato et al. use the uniform 2D GP 
equation for a small u± limit: 

"dT ~ \ 2M 

To confirm whether the cascade is direct or inverse, the thermalization pro- 
cess of decaying turbulence in an isolated system is an effective indicator. 
If the cascade is direct, an essential part of the energy reaches the small- 
est scales available in the simulation and the system quickly evolves toward 
thermodynamic equilibrium filled with short scale excitations of fluid. On 
the other hand, if the cascade is inverse, fluid motions of the system size are 
strongly excited even at later stages. To check this thermalization process, 
Numasato et al. introduce neither dissipation nor energy injection in this 
work. 

The method to generate turbulence is almost the same as that used for 



3D decaying turbulence [68] : the initial condition of \P is set to constant 
density p± and random phase 0, which varies at large scales. From this 
initial condition, we obtain, during time evolution, a 2D QT composed of a 
random configuration of quantized vortices. A simulation was performed in 
a periodic box with a size of 64£ for all the results shown below. 

Figure HH shows the evolution of the total energy E, the kinetic energy 
-Ekin, and its compressible and incompressible parts E^ in and -E kin , and the 
number of vortices. E, E^ n , E^ in , and E l kin are obtained from Eqs. f|4"T|) . 
fl4"8"j) . and fpEl by replacing \1/ and p with ty± and p±. Because the system 
has no dissipation, E is time independent. At the initial stage (t < 5), an 
intensive process of vortex creation leads to fast transformation of E£ in into 
El in , which is larger than E^ in for the time interval 2 < t < 4. At later 
stages (t > 5), E^ in is practically independent of time. The largest maximal 
value of the vortex number is achieved at the crossover time t ~ 5. Its decay 
is faster for larger g± (not shown), because the probability of the dominant 
nonlinear process of vortex pair annihilation is larger for larger g±. A decay 
of the vortex number leads to an increase in the flow compressibility. The 
compressibility measured by E^ in /El in is larger for larger g± (not shown). 
The system finally reaches its equilibrium state with no quantized vortex. It 
takes longer to reach its thermodynamic equilibrium state for smaller g±. 
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Figure 14: Time evolution of E, Skim -^kin' ^kiw an< ^ the number of vortices. A simu- 
lation was performed with gj^ = 4, and with length, energy, and time normalized by 
H 2 /(2M£j_), and (M£j_)/h, respectively. [Numasato, Tsubota and L'vov: Phys. Rev. A 
81 (2010) 063630, reproduced with permission. Copyright 2010 the American Physical 
Society] 



The energy spectra E^ n (k) and E£ in (k) at different moments of time are 
shown in Fig. [151 Here, E^ n (k) and E^ in (k) are defined by 

«W=2(2^/ d ^(« 1 ) 2 ' < 86 ' 

for [pjj ' 1 = [a/pT^s] ' 1 and its Fourier transformation p c _£, and the 2D total 
number of particles N± = J dr\^±\ 2 . At around t ~ 3, E l kia {k) are close to 
the Kolmogorov law E kin (k) oc A; -5 / 3 as shown in Fig. [T5] (a). This behavior 
is related to the energy cascade and will be discussed later. At later times, 
the energy begins to accumulate at large k and the E^ n (k) asymptotically 
approach the quasistationary state. E^ in (k) at this stage are close to the 
thermodynamic equilibrium with energy equipartition between degrees of 
freedom. In 2D systems, this gives E kin (k) oc A; as shown in Fig. [T51 (c). This 
will be also discussed later. This distribution, however, does not correspond 
to that under thermodynamic equilibrium. This is related to the fact that 
the system does not achieve full equilibrium at these times. Indeed, Fig. [TH 
(c) shows that E^ n continues to converge to E kin . We interpret this stage as 
a kind of flux equilibrium, when the E kin (k) is determined by the energy flux 
from El in to E£. in . The exchange between E\.- m and E kin plays a subdominant 
role at these times. After a long time evolution, most of E^ n finally comes 
to comprise E kin . This is the full thermodynamic state. As expected, E kin (k) 
is proportional to k as shown in Fig. [15] (c). 

Numasato et al. next consider the particle number spectrum N(k) = 
|\I/_i_| 2 = pj_, where is the Fourier transformation of Figure [TBI (a) 
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Figure 15: (a) E^ in (k) at an earlier moment of time, (b) E^^k) a ^ a later moment of 
time, (c) E^ in (k) at a later moment of time. In all figures, g± — 4. [Numasato, Tsubota 
and L'vov: Phys. Rev. A 81 (2010) 063630, reproduced with permission. Copyright 2010 
the American Physical Society] 



shows N(k) for different time intervals. Through the turbulent state, the 
spectrum asymptotically approaches the form N(k) oc k^ 1 . To rationalize 
this behavior, we note that this dependence should be related to E^ m {k) oc k 
in the full thermodynamic equilibrium. Indeed, in this state, the interaction 
energy can be neglected in comparison with the kinetic and quantum energy. 
In addition, the quantum energy corresponds to zero-point motion and does 
not give rise to particle currents. Thus, only the kinetic energy spectrum is 
related to the particle spectrum. Moreover, almost all E^in becomes E^ in and 
there are no vortices, and as a result, has no singularity and becomes of 
order unity. Thus the kinetic energy density can be written as follows: 



= f drp ± m* ~ — -g- / dk k^\ 



2(2n) 2 N ± M 2 
As a result 



2N ± M 2 J r 1 ri 2(2tt) 2 N ± M 2 

h 2 , 

dk k 3 p±. 



(87) 



E< (k) ~ h2k2N{k) (88) 
hkin[k} ~ 2(2k) 2 N ± M 2 - m 

In this way, the two relations E^ in (k) oc k and N(k) oc k~ x hold consistently. 
This relation is similar to the relation between enstrophy Q and kinetic energy 
E in 2D CT. 
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Figure 16: (a) Particle number spectra, (b)-(c) Averaged incompressible kinetic energy 
flux in a short time period r = 0.20 in the time range 2.50 < t < 4.95. In all figures, 
gj_ = 4. [Numasato, Tsubota and L'vov: Phys. Rev. A 81 (2010) 063630, reproduced 
with permission. Copyright 2010 the American Physical Society.] 

Next, Numasato et al. consider the value and direction of the compress- 
ible and incompressible energy flux IT£ in and rP kin . Numerical results for II£ in 
and n^jj are shown in Fig. [16] (b) and (c). Ilj^ takes positive values for 
3(2tt/L) < k < 2vr/£ ± at least for 2.50 < t < 4.95. This strongly supports 
the idea that ££ in propagates from small k to large k and we conclude that 
at some region of the system parameters Numasato et al. can observe a 
2D-direct energy cascade. 

When the Kolmogorov spectrum is formed, a 2D Richardson cascade can 
be seen. Numasato et al. choose the shortest intervortex pair length Z p for 
all vortices. In Fig. [17] (a), the averaged vortex pair number is shown and 
is proportional to l~ n , where n depends on g± and 1.30 < n < 2.12. This 
power law suggests a self-similar spatial structure. For 2D vortices, one of 
the most effective lengths, corresponding to the length of 3D vortex ring, is 
the intervortex length. 

Important information about the motion can be extracted from the fre- 
quency power spectrum, which is the Fourier transformation of the different- 
time pair correlation function. Numasato et al. observe this frequency power 
spectrum for the compressible velocity component: 

J dte^vt in (k,t').{vl in (k,t' + t)y, (89) 

at later times. As the system evolves to the thermodynamic equilibrium 
state, the power spectrum forms sharp peaks. In Fig. [17] (b), Numasato 
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Figure 17: (a) Averaged vortex pair numbers as a function of intervortex length l p . (b) 
Position of the maximum frequency power spectra of the compressible velocity component 
for different wave vectors, averaged over a long time in the state of full thermodynamic 
equilibrium (full dots). Bogoliubov's frequency spectrum w max , Eq. (|90)) . and its large-fc 
asymptotic form of w max oc k 2 , where the wavenumber, energy, and frequency are normal- 
ized by h 2 /(2M£l_), and (M^J/tL In both figures, g ± = 4. [Numasato, Tsubota 
and L'vov: Phys. Rev. A 81 (2010) 063630, reproduced with permission. Copyright 2010 
the American Physical Society] 



et al. plot the position of the maxima integrated over a long time interval 
for t in the thermodynamic equilibrium state. The eigenfrequencies of the 
thermal fluctuations have been determined by Bogoliubov to be 



/ h 2 k 4 g±P±k 2 

"w ^Viii^ + ^w - (90) 

The excellent agreement between the theoretical and numerical results indi- 
cates that the observed thermal fluctuations of the compressible velocity com- 
ponent do indeed correspond to Bogoliubov's elementary excitations. The 
relatively small but finite width of the peak characterizes the finiteness of 
the lifetimes of these fluctuations, caused by interaction of the fluctuations 
with different k. 

Numasato et al. finally note that E^ in approaches and there is no vortex 
in the thermodynamic equilibrium state. The temperature of this state is, 
therefore, below the Kosterlitz-Thouless transition temperature Tkt- If we 
increase the initial energy injection, we expect the final equilibrium state with 
many randomly paired nucl eated and annihilated vortices, the temperature 



of which is above Tkt |l08j |. We also note that 2D QT is not restricted to 
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theory. For atomic BECs, the initial condition can be prepared by the phase 
imprinting technique. By using experimentally realistic parameters M = 
1.46 x 10" 25 kg, a = 5.77nm, N = 10 3 , a z = 1.34/im, a± = y/h/(Mu ± ) = 
4.25/im, the healing length is estimated to be £j_ = 1.22//m. As a result, the 
size of the system in x-y space must be L ~ 2a_|_ ~ 7£ to see the effects 
discussed in this section. 



3.5.5. Quantum turbulence in atomic Bose-Einstein condensates 

The study of the turbulent state in quantum fluids and its relation to 
CT is an intriguing physical problem. Although the study of QT has a long 
history, only superfluid 4 He and 3 He systems have been used to realize QT 
until recently. Recently, atomic BECs have become another can didate for 
QT resea rch, since a turbulent state was realized in this system [lO^, |n3, 



111, 112 



Compared with a helium system, the characteristics of trapped BECs are: 
(i) a BEC system is weakly interacting and can be easily treated theoreti- 
cally, (ii) many physical parameters of BECs are experimentally controllable, 
and (iii) various physical quantities such as the density and phase of BECs 
can be directly observed. Quantized vortices can be considered to be holes 
of density and singularities of phase. Shortly after trapped BECs were first 
realized, experimental groups reported vortex lattice structures, as well a s 



the crystallization dynamics of these structures under rotation [1131 . Il 14 
These dynamics have been successfully con firme d quantitatively by numer- 
ical simulations using the GP equation 70l. Il 151 ] . However, in experimental 



research on trapped BECs, another important phenomenon of quantized vor- 
tices, namely QT, has not been adequately studied until recently. Noting that 
quantized vortices are observable and that almost all physical parameters of 
trapped BECs are controllable, such systems are an ideal prototype for truly 
controllable QT. QT in trapped BECs is, therefore, used to determine several 
details of the system, such as the distribution of vortex length, details on the 
cascade of vortices, the isotropy or anisotropy of vortex configuration, and 
details on correlations among vortices related to eddy viscosity, as already 
considered for CT [i| . Clarifying any of these will lead to the detailed under- 
standings of the transition to QT and its universality. Therefore, research 
into QT offers the promise of greater advances in understanding turbulence 
than has been possible in past studies of turbulence. 

There is a disadvantage in using trapped BECs to study QT: to generate 
turbulence, we cannot apply a velocity field, which is widely used for research 
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on CT and QT of superfluid helium, because BECs are trapped. To realize 
turbulence in this system, several ways have been theoretically proposed 
such as relaxation from a s tron gly degenerate nonequilibrium gas across the 
BEC critical temperature jl!6l | or crystallization of an i solated BEC from 



the vortex-free to vortex lattice state under rotation jl!7 



Here, we present a precession rota tion which has two rotation axes, one of 



which rotates around the other [118| . We start from the GP equation under 
the rotating field fl: 

(i - 7) h— = — + g |*| 2 - n + V - ihQ ■ r x V J *. (91) 



dt V 2M 
Here V is the trapping potential satisfying 

V = K 1 - W ~ SyW + (1 + W - S y )y 2 + (1 + 5 y )z 2 } , (92) 

with the trapping frequency u and elliptical deformation parameters 5 Z and 
5 y in the x-y and z-x planes. To develop the BEC to a turbulent state rather 
than a vortex lattice state, we use precession rotation along the z and x axes; 
the first rotation along the z axis rotates around the second rotation along the 
x axis. The resulting rotation field becomes fl = (Q x , Q z smQ x t, fl z cos fl x t), 
where fl z and fl x are the frequencies of the first and second rotation, respec- 
tively. The advantage of using this precession rotation to study turbulence 
is high controllability of the state from a nonturbulent vortex lattice to fully 
developed turbulence by changing the ratio Q x /fl z . 

Now we consider a system at very low temperatures. To apply the re- 
sult shown in Fig. [13] to the dissipation term 7, we consider the Fourier 
transformed form of Eq. 



(* - 7)ft-^r = —77 - a* ) * + 9 Y + v - inn ■ R, 



dt \2M 



(93) 

R= I dr e~ ikr r x V*, 



where 7, Y, and V are defined in Eqs. f[53|) and (jog)) . In this work, we sub- 
stitute the result at T = 0.01 shown in Fig. [13]into 7 with the healing length 
at the trap center: £ = h/ \j2Mgp{r = 0) at t — 0. For other numerical 
parameters, we use the following, taken from experiments on 87 Rb atoms: 
M = 1.46 x 10" 25 kg, a = 5.61 nm, N = 2.50 x 10 5 , and u = 150 x 2tt Hz. 
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We start from a stationary solution without rotation and elliptical defor- 
mation. At t = 0, we turn on the rotation Q x = Q z = 0.6u and elliptical 
deformation 5 Z = 5 y = 0.025, and numerically calculate the time develop- 
ment of the GP equation fl93|) . In the initial stage, vortices start to enter the 
BEC making the system quite anisotropic. After tu ~ 150, the BEC recovers 
isotropy and the system enters a statistically steady state. The steady turbu- 
lence is sustained by the balance between the large-scale energy injection due 
to the rotation and the small-scale dissipation. Furthermore, in all stages of 
the dynamics, E kin is always much larger than E^ in and the dynamics of the 
BEC are dominated by vortices rather than compressible excitations. 

Figure [TS] (a) shows E l khi {k) and ITj^ for the steady turbulent state. 
E kin (k) satisfies the Kolmogorov law in the inertial range 2n / R^-p < k < 
2?t/£, where Rtf = \f2p{t = 0)/(Mu) is the Thomas-Fermi radius and 
represents the largest scale in the BEC. Furthermore, the energy flux is 
nearly constant IT^ ~ 1 Afou 2 / '(N M) in the inertial range, supporting the 
fact that the incompressible kinetic energy steadily flows in wavenumber 
space through the Richardson cascade at the constant energy transportation 
rate ej^ = rT kin . Using this e j^ , we can estimate the Kolmogorov constant 
C ~ 0.25, which is smaller than that in CT and consistent with the work for 
the uniform system in Sec. 13.5.21 

To investigate the relation between the Kolmogorov law and the Richard- 
son cascade, we calculate the vortex length distribution n(l)Al inside the 
condensate, where n(l) Al represents the number of vortices with length from 
/ to / + Al. As shown in Fig. [18] (b), n(l)Al obeys the scaling property 
n(7) oc l~ a for 27r£ < / < 27ri?TF- This reflects the self-similar Richar dson 
cascade in which large vortices entering the condensate from the surface 115 
are divided into smaller vortices. The scaling exponent a is close to unity, 
which is consistent w ith t hose given by Araki et al. (a ~ 1.34) [83] and 
Mitani et al. (a ~ 1) [ll9l 



To visualize the turbulence, we plot the isosurface of p and the spatial 
distribution of the vortices inside the condensate in Figs [19] (a)-(f). At 
tu = 10, the surface of the BEC becomes unstable (Figs. [[9] (a) and (d)), 
and vortices appear in the BEC at tu = 50 (Figs. [[9](b) and (e)). Figures fT9l 
(c) and (f ) shows QT with no crystallization but with highly tangled vortices 
at tu = 300. 

Finally, we note that the obtained energy spectrum in Fig. [18] (a) is not so 
clear straight line and its consistency with the Kolmogorov law is incomplete. 
This inconsistency comes from the anisotropy of turbulence around the y- 



52 



10° 
lO- 2 

J io- 4 

10" 6 



io- s 



0.1 



(a) 

2tt/% f 2tt/£ 



• • 




> 




Cfe 2 / 3 fc" 5 / 3 \ 

(C=0.25) * - 




"kin 1 

1 i ■ 



10 



1.5 



100 



1 10 



o.i 



(b) 

2tt£ 2tt% f 



A;=2.5 




Figure 18: (a) £^ in and II^ in . (b) Vortex length distribution n(l)Al inside the Thomas- 
Fermi radius Rtf- In both figures, length, wavenumber, energy, and time are normalized 
by ah = y/H/(Mv), hw, and l/v respectively. E^ in , 11^, and n(l)Al are obtained 

from an ensemble average of 25 randomly selected states at t > 300. [Kobayashi and 
Tsubota: Phys. Rev. A 76 (2007) 045603, reproduced with permission. Copyright 2007 
the American Physical Society] 



axis around which there is no rotation, and is improved by other simulations 
of Q T of trapped BECs under the precessional rotation around three axes 

fiia j. 

Recently, a turbulent state has been realized in atomic BECs using two 
methods. Weiler et al. perfor med a rapid quench of an 87 Rb gas through the 
BEC transition temperat ure [109j| , which is similar to the method advocated 



by Berloff and Svistunov [121| . Through the high density fluctuation regime 
(weak turbulence) in a short period, several vortices and anti-vortices were 
formed, creating the system turbulence. As a method with better control 
of the turbulence, the experimental group o f Begnato intr oduced an exter- 



nal oscillatory perturbation to a 87 Rb BEC llOl . IlllL Ill2j ]. This oscillating 



magnetic field was produced by a pair of anti-Helmholz coils which were not 
perfectly aligned to the vertical axis of the cigar-shaped condensate. Addi- 
tionally, the components along the two equal directions that result in the 
radial symmetry of the trap were slightly different. This oscillatory field 
induced a coherent mode excitation in a BEC. For small amplitudes of the 
oscillating field and short excitation periods, dipolar modes, quadrupolar 
modes, and scissor modes of the BEC were observed, but no vortices ap- 
peared. Increasing both parameters, the vortices grew in number, eventually 
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Figure 19: (a)-(c) Isosurface plots ol 5% of the maximum condensate density p. (d)- 
(f) Configuration of quantized vortices inside the Thomas-Fermi radius Rtf- (a), (d) 
too = 10, (b), (e) tu> = 50, (c), (f) tcu = 300. The method for identifying vortices in (d)-(f) 
is the same as that in Fig. M [Kobayashi and Tsubota: Phys. Rev. A 76 (2007) 045603, 
reproduced with permission. Copyright 2007 the American Physical Society] 



leading to the turbulent state. In the turbulent regime, they observed a rapid 
increase in the number of vortices followed by proliferation of vortex lines 
in all directions, where many vortices with no preferred orientation formed 
a vortex tangle. Another remarkable feature is that a completely different 
hydrodynamic regime followed: the suppression of aspect ratio inversion dur- 
ing free expansion, despite the asymmetric expansion (from a cigar-shaped 
to pancake-shaped) of the usual quantum gas of bosons, or the isotropic ex- 
pansion of a thermal cloud. Although the theoretical understanding of this 
effect remains incomplete, it represents a remarkable new effect in atomic 
superffuids. 

For a better understanding of the experiment al re sults, we performed 



numerical simulations based on the GP equation [112j ]. The net potential 



V acting on the atoms is the sum of the harmonic magnetic trap and the 
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oscillatory field, and can be approximately expressed by 

M 2 
V — — [u x {x cos^i + y sin^i — z sin# 2 — 6i(l — cos fi t)} 

+ u 2 {ycos6 1 -xsin^ - 6 2 (1 -cosQ t)} 2 ( 94 ) 
+ u 2 {z cos 8 2 + x sin 6 2 — 6 3 (1 — cos Qot)} 2 ] , 

where 9i = — cos^o^) are time dependent angles. For the experimental 
conditions, iV = 3 x 10 5 , w x = 2vr x 23Hz, w r = 2vr x 210Hz, fi = 2vr x 200Hz, 
Ax ~ 7r/60, and A 2 — 7r/120 were used. The amplitudes for the translational 
oscillation of the potential minimum are (61,62,63) = aa r (2, 5, 3)/im, where 
a r = y/h/ (Moj r ) and a is a variable parameter that represents the ampli- 
tude of the center-of-mass oscillation, being proportional to the amplitude 
of the excitation. For simplicity, we employ \I> = *& r (y, z)^f x (x) and consider 
2D simulations in y—z space. Here, we consider the BEC surrounded by a 
thermal cloud and use the constant 7 for the GP equation f[3"2"j) . Since the 
thermal atoms, which are the origin of the dissipation, move together with 
the potential, we also have to consider the reference frame co-moving with 
the potential. In this frame, the GP equation becomes 



(i-nf)h 



ft 2 V r 2 
' 2M 



+ V r -/1 + g r \^r\ 2 ~ fi(t) • L - V(t) ■ P 



*r, (95) 



with momentum P = —iKV and angular momentum L = —ihr x V. Here 
V 2 = d 2 /dy 2 + d 2 /dz 2 and g r = Akcl/R x , with R x being the characteristic 
size of the condensate along x axis. We consider L = (fl x , 0, 0) sinf2 £ and 
v = (0, v y , 0) sin Qot- Using half of the oscillation period T = tt/Qq, we 
obtain v y ~ 2<5 2 /T = 2Q 62/tt. The rotation frequency Q x is also estimated 
as fl x ~ 2y4 2 /T = Q /Q0, providing a very small contribution. 

Figure |20] shows snapshots of the density profile for different excitation 
times ranging from 13 to 17ms. Additionally, we have calculated the mean 
angular momentum per atom, (L x ) = J dr ^>* r L x ^> r , as a function of the ex- 
citation time. Using a = 1.6 and 7 = 0.02, our simulation shows that (L x ) 
blows up after 15 ms of excitation. At this point, the condensate forms wavy 
patterns which develop to dark solitary wave s wh ich subsequently decay into 



several vortex pairs via the snake instability |122| . A more complex dynamic 
takes place after the first events of vortex formation, consisting of the gener- 
ation of an undetermined number of vortices, characterizing the emergence 
of the turbulent regime. As a increases, the nucleation of vortices occurs 
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Figure 20: Snapshots of the BEC after different times of excitation [112| . Figure shows 
the 2D plot of the density profile. The colors range from red (high density) to blue (low 
density). [Seman, Henn, Shiozaki, Roati, Poveda-Cuevas, Magalhaes, Yukalov, Tsubota, 
Kobayashi, Kasamatsu, Bagnato: Laser Phys. Lett. 8 (2011) 691, reproduced with per- 
mission. Copyright 2011 WILEY- VCH Verlag GmbH k Co. KGaA.] 



at earlier times, with a faster evolution to QT. This agrees well with the 
observations. Furthermore, the time scale of the vortex events in the simula- 
tion, which is of the order of 10 ms, is consistent with the times observed in 
the experiment. These results demonstrate that the combination of rotation 
and translation is essential to produce vortices. The dynamics of the BEC 
strongly depend on the strength of the dissipation 7. If the dissipation is 
absent, no instability associated with soliton creation occurs. Values of 7 
between 0.015 and 0.025 are optimal for the generation of vortices and QT. 
The simulations cannot reproduce the full experimental results, since the 
experimental system is a 3D gas. Nevertheless, good qualitative agreement 
with the experiment has been achieved. 

3. 6. Cascade process in quantum turbulence 

The most important concept for understanding QT is the cascade process 
of the energy and vortices as well as CT. In Q T, there are two regions of 



the cascade process in wavenumber space [123J. The first region is called 
the classical region below the inverse of the mean intervortex spacing. The 
dynamics of vortices in the classical region are dominated by the Richardson 
cascade, in which large vortices are broken up self-similarly into smaller ones, 
or the collective dynamics of aggregated quantized vortices at scales larger 
than the intervortex spacing. Such behavior of vortices supports the analogy 
of QT to CT, namely the Kolmogorov energy spectrum [5j. The second region 
is called the quantum region, in which vortex dynamics are dominated by the 
effects of the quantized circulation, sp ecifically the Kelvin wave cascade of 



vortices, which does not appear in CT 124 . Il25| . The Kelvin wave cascade 
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is also a very important concept in understanding the dissipation mechanism 
of QT at very low temperatures. 

Here, we briefly summarize the overall picture of the energy spectrum of 
QT at zero temperature pj based on theoretical and numerical studies (Fig. 
2~Tj) . If a vortex tangle in QT is homogeneous and isotropic, there are two 
characteristic length scales: the mean intervortex spacing I = L~ 1//2 with a 
vortex line length density L, and the healing length £ corresponding to the 
size of the vortex core. I is much larger than £ in the case of superfluid 
helium, while both are usually of the same order for atomic BEC. Here, we 
consider the former case for simplicity. Using / and £, we can define the 
corresponding wavenumbers k\ = 2n/l and = 27r/£. At length scales 
larger than /, the dynamics of QT are dominated by a tangled structure of 
many vortices. Because vortex dynamics become collective at large scales, 
quantization of the circulation is not relevant and the dynamics are similar to 
those of eddies in CT. This is why this region can be referred to as the classical 
region. As a result, the energy spectrum E(k) in the range k < k\ obeys 
the Kolmogorov law. In the classical region, vortices sustain a Richardson 
cascade that transfers energy from smaller to larger wavenumbers without 
dissipation. The Richardson cascade can be understood as large vortices 
breaking up into smaller ones in real space. 

Vortices in QT can reconnect many times (see Fig. which is the domi- 
nant dynamics at length scales comparable to I. Through the reconnections, 
small cusps or distortion waves are formed on the vortex l ines , whi ch are 
regarded as the primary source of Kelvin waves in QT 32|, 1 124 . Il26| . The 



wavelength of the created Kelvin waves is of the order of I. 

At length scales smaller than I, which is referred to as the quantum region, 
the Richardson cascade is no longer dominant and the quantize d circulation of 
vortices and motion of each vortex line become significant QjJ, 1 124 . Il27l . Il28 . 
129l . Il30l . Il3ll . Il32l . Il33l . Il34j | . In this range, vortex dynamics are characterized 
by the cascade process of the Kelvin waves formed by reconnection. The 
nonlinear interaction of the Kelvin waves is the origin of the cascade from 
small to large wavenumbers. The energy spectrum in the quantum region 
ki < k < k^ is theoretically predicted to obey a Kolmogorov-like power law: 
E(k) oc k v . Three th eoretical va lues of 77 have been pr edicted: 7 7 = —7/ 5 



^,[3 SEES, V = - 5 / 3 HO, and 77 = -1 [l26l . [Tiil . fl3ll . [l32 



At finite temperatures where the mutual friction between superfluid and 
normal fluid is effective, Kelvin waves and relevant turbulent flow are strongly 
dissipated by the viscosity of the normal fluid, and the cascade of Kelvin 
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logE(fc) 




Figure 21: Overall picture of the energy spectrum of QT at zero temperature. The energy 
spectrum depends on the scale and its properties change at about the scale of the mean 
intervortex spacing I. When k < ki = 2n/l, a Richardson cascade of quantized vortices 
transfers energy from large to small scales, maintaining the Kolmogorov spectrum E(k) = 
Ce 2 / 3 fc~ 5 / 3 . When k > ki, energy is transferred by the Kelvin-wave cascade, which is 
a nonlinear interaction between Kelvin waves of different wavenumbers. In this region, 
the energy spectrum also takes the power-law structure E(k) oc k~ v and the three values 
T) = —7/5, rj = —5/3, and r) = —1 have been predicted. Eventually, energy is dissipated 
at scales of £ by the radiation of elementary excitations. The crossover region between 
Richardson cascading and Kelvin-wave cascading regions remains an open question, and 
two scenarios, E(k) oc k 2 and E(k) oc /c _3 -fc°-fc _1 , are proposed (see text). 



waves turns off. 

In the region k ~ k^, 
elementary excitations, such as phonons and rotons, in the primary decay 



ke, the Kelvin waves with wavelength £ change to 



process of QT near zero temperature [125 



There is one open question regarding the energy spectrum in the region 
k ~ ki, namely the transitional region between the Richardson cascade and 
the Kelvin-wave cascade, referred to as the classical-quantum crossover. A 
theoretical study proposed a bottleneck effect connecting the spectrum sat- 



isfying the thermalization spectrum E(k) oc k [1361 ] . Another theoretical 
prediction is based on vortex reconnection dynamics on the scale ~ I, 



in 



which the transitional casc ade process was predicted to occur by reconnec- 



tion of the vortex bundle 137|. The crossover range is divided into three 
subranges, giving E(k) oc k~ 3 , E(k) oc k°, and E(k) oc k^ 1 in each region. 
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We discuss this spectrum in Sec. 13.6.31 



3.6.1. Classical region 

There is a significant issue regarding the energy spectrum in the classical 
region k < k\ in terms of the analogy of QT to CT. In this region, the energy 
spectrum is determined by the collective behavior of many vortices, such as 
the vortex tangle and the aggregated bundle structure at scales larger than 
/. Several numerical studies have calculated the energy spectrum in this 
region by simulating QT at zero temperature. We have numerically found 
and already discussed the Kolmogorov energy spectrum (146]) by using the 
GP model in Sec. 13.51 (Sec. I3.5.2[) . Araki et al. also found the Kolmogorov 



energy spectrum through the vortex-filament model [83 



In our numerical studies discussed in Sec. 13.5.2] however, the system size 
was not so large and the inertial range was less than one order in wavenum- 
ber space. Furthermore, the mean intervortex spacing I was close to the 
healing length £, being too short to study the Kelvin wave cascade. To ob- 
tain the energy spectrum of a wider range of wavenumber space, Yepez et 
al. performed a large-scale simulation of th e GP model by using a novel 



unitary quantum lattice gas algorithm [1381 1 . They found that the incom- 



pressible kinetic spectrum E^ in had three distinct power-law k~ a regions 
that ranged from the classical turbulent regime of Kolmogorov a = 5/3 
at large scales k < (a/3/(27t))L/^ to the quantum Kelvin-wave cascades 
a = 3 at small scales k > (v / 3/2)L/^. There was a semiclassical re- 
gion 6.34 < k < 7.11 connecting the Kolmogorov and Kelvin-wave spectra 
(V3/(2n))L/£ < (y/3/2)L/£. Compared with our simulations, this simula- 
tion supplied the Kolmogorov spectrum over a much wider inertial range, 
with about two orders in wavenumber space. Although they related the 
k~ 3 spectrum at small scales k > (\/3/2)L/£ to the Kelvin-wave cascade, 
the length scale in this region is smaller than the vortex core size. Hence 
the k~ 3 spectrum comes most probably not from the Kelvin-wave cascade 
but f rom the velocity profile abour a vortex, as pointed by several authors 



1391. 1141 1141 



To clarify the Kolmogorov energy spectrum in the classical region and 
the energy spectrum in the classical-quantum crossover, Sasa et al. also 
performed a simulation of the GP equation at much larg er sc ales than that 



discussed in Sec. 13.5.21 i.e., from L = 128£ to 512£ [1421 ] . In contrast 
to the work by Yepez et al, Sasa et al. focused on the classical region 
(and classical-quantum crossover)of k < I . We numerically solved the GP 
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Figure 22: Left: Simulation results for the incompressible kinetic energy spectra E^ in (k). 
A varies from A ~ 1.5 for L = 128 to A ~ 2.2 for L = 512. Dot-dashed line: Kolmogorov 
spectrum. In the figure, length and energy are normalized by £ and fi 2 /(4M( 2 ). Right: A 
snapshot of vortex lines at the fully developed turbulent state of L = 512 demonstrating 
the self-similarity of the bundle structure (see the dotted circles representing the zoomed 
regions whose vortex distributions are shown subsequently), typical for fully developed 
turbulence. [Sasa, Kano, Machida, L'vov, Rudcnko, and Tsubota: Phys. Rev. B 84 (2011) 
054525, reproduced with permission. Copyright 2011 the American Physical Society] 

equation (|52|) with the dissipation term 7 of the step function form (p>4"l) with 
7o = 1. We investigated the decaying turbulence without energy injection, 
and a uniform density p — 1 and a random as the initial wave function 
to create our system turbulence. The random phase is arranged in the same 
way as in our previous work discussed in Sec. 13.5.21 for decaying turbulence, 
i.e., placing random numbers between — an to an at every distance A and 
connecting them smoothly, where a is the control parameter for the energy 
injection. 

The main numerical results are shown in Fig. [22j The left panel shows the 
incompressible kinetic energy spectrum E^ in (k). The Kolmogorov spectrum 
extends to the lower-fc range and increases with the system size. The visible 
extent of the Kolmogorov spectrum is much larger than that in all the results 
including our previous simulations. The right panel of Fig. [22] displays large 
self-similar structures of tangled vortices in the fully turbulent state: large- 
scale vortex bundles in the maximum size, 512£, and smaller self-similar 
tangled structures inside this cubic region in the subsequent insets. The 
visualization of vortices clearly shows the bundle structure, which has never 
been confirmed in GP simulations in smaller boxes. 

An important observation is a plateau-like region for k > 1.5£, which is 
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a definite pileup over the Kolmogorov spectrum; this is a clear manifesta- 
tion of energy stagnation. We exp ect t his plateau-like structure to be an 
indication of the bottleneck effect 136| proposed by L'vov et al. in the 
classical-quantum crossover region, which will be discussed in Sec. 13.6.31 



Detailed comparisons are discussed in Ref. |142 



3.6.2. Quantum region 

In the region of k > ki, the picture of aggregated vortices is no longer 
effective and the motion of each vortex line becomes essential. The most 
probable dominant dynamics of vortices are Kelvin waves, which originate 
from distortion waves on the vortex lines after their reconnection. A Kelvin 
wave is a transverse, circularly polarized wave motion, with the approximate 
dispersion relation for a rectilinear vortex: 



nk 2 
Arc 



+ c 



(96) 



with a dimensionless constant c ~ 1. A; is the wavenumber of the Kelvin 
wave, and is different from that used for the energy spectrum. Kelvin waves 
were first obse rved by inducing torsional oscillations in a rotating superfluid 
4 He [l43l.[l44|. 

Although the Kelvin-wave cascade seems to be a very important mech- 
anism in QT at scales smaller than I at very low temperatures, it is a non- 
trivial problem regarding the actual cascade process. At finite temperatures 
where there is a significant fraction of normal fluid, Kelvin waves are damped 
by mutual friction. On the other hand, at very low temperatures they can 
be damped only by the radiation of phonons. Vinen estimated the rate of 
radia tion, and found that it is extremely low unless the frequency is very 
high 125|. In QT, the main origin of the Kelvin- wave nucleation is the vor- 



tex reconnection. When two vortices reconnect, they twist to become locally 
antiparallel at the reconnection point and create small cusps or kinks after 



the reconnection which were confirmed numerically [28l . l29l |32| . Svistunov 



suggested that the relaxation process of these cusps or kinks causes the emis- 
sion of Kelvin waves, and plays an important role in the decay of QT at low 
temperatures. Following Svistunov, Vinen et al. performed a numerical sim- 
ulation of a Kelvin wave excited alon g a si ngle vortex line using the vortex 
filament model discussed in Sec. 13.21 1281 1 . The authors consider a model 



system in which helium is contained in the space between two parallel sheets, 
separated by a distance £b = 1cm, with a single, initially rectilinear, vortex 
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stretched between opposite points on the two sheets. Kelvin waves can be 
excited on this vortex, and periodic boundary conditions are applied at each 
end. The allowed wavenumbers of the Kelvin waves are given by 



k 



2nn 



(97) 



where n is a natural number. The authors imagine that the mode with a small 
integer n is continuously driven. As the amplitude increases nonlinearly, 
coupling to other modes occurs and they can expect energy to flow from 
the mode no to other modes with both large and small wave numbers. To 
mimic the effect of phonon emission, they introduce strong damping for all 
modes with n exceeding a larger critical value n c . Then, they can obtain 
a statistically steady state in which the energy injection in the mode with 
no is balanced by dissipation in the modes with n > n c and calculate the 
corresponding energy spectrum. The authors observe no reconnection. 

The simulations are based on the vortex filament model with the full 
Biot-Savart law based on Eq. ( II 2p . The force that drives one mode is of 
the form Vpnsm(k z — co t), where k = 2nn /£B, p is the density of the 
helium, and Uq is related to ko by the dispersion relation (|96|) . Damping 
at the highest wavenumber k c = 27rn c /-^B = l/60cm _1 is applied using a 
periodic smoothing process. We calculate the root mean square amplitudes 
Ck(t) = (Ck(k) 1 ^ 2 of the Fourier components of the displacement of the vortex 
from the original straight line. Figure [23] (a) shows how these amplitudes 
develop in time after the application of a drive with V = 2.5 x 10 _5 cm s _1 
and ko = lOvrcm^ 1 . We see that initially only the mode of ko is excited. 
However, as time passes, nonlinear interactions lead to excitation of all other 
modes. Eventually the spectrum reaches a statistically steady state. For 
large values of k, where the modes practically form a continuum, the steady 
state is observed to have, to a good approximation, a spectrum of the simple 
form 



where the dimensionless parameter A is of order unity. 

Figures [23] (b) and (c) show the effects, respectively, of increasing the 
drive amplitude V by a factor of 10 and of changing the drive wavenumber 
ko- We see that there is little effect on the steady state, within the error 
of the simulations, at least at the higher wavenumbers. The steady state 



Ck — Ai B x k 3 



(98) 
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(a) (b) (c) 
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fc(l/cm) *(l/cm) fe(l/cm) 

Figure 23: (a) Time development of Ck(t) in the Kelvin- wave cascade. The short-dashed, 
dash-dotted, dotted, long-dashed, and solid lines refer, respectively, to averages over 0-800, 
10000-10800, 20000-20800, 40000-40800, and 140000-140800 s. (b) Steady state values of 
£fc for two different drive amplitudes. The solid line and dotted line are for, respectively, 
V = 2.5 x 10~ 5 cm s _1 and V = 2.5 x 10~ 4 cm s -1 . The long-dashed line has the form 
of Eq. (|98p. (c) Steady state values of for three different wavenumbers. The dotted, 
short-dashed, and solid lines refer, respectively, to ko = 2-7T cm , fco = 47rc m , and 
fco = 107r cm -1 . Again the long-dashed line has the form of Eq. (|98|) . [Vinen, Tsubota, 
and Mitani: Phys. Rev. Lett. 91 (2003) 135301, reproduced with permission. Copyright 
2003 the American Physical Society] 



takes longer to be established at the lower drive amplitude, which suggests 
that even with a small drive amplitude, the same steady state would be 
established after a sufficiently large time. 

The mean energy per unit length of a vortex in mode k is related to (k 
by the equation 

E K (k) = e K k 2 C 2 k , (99) 
where ex is an effective energy per unit length of vortex, given by 




(100) 



It follows from Eqs. ® and (|22} that 

E K (k) = Ae K (k£ B )-\ (101) 

Thus the steady state is characterized by the energy spectrum (llOip . and 
this spectrum is insensitive to the frequency and amplitude of the drive and 
to the power input at the drive frequency. 
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What happens to this kind of Kelvin-wave cascade in QT? When two 
vortices reconnect, they twist to become locally antiparallel at the reconnec- 
tion point and create small cusps or kinks after the reconnection. Svistunov 
suggested that the relaxation of these cusps or kinks causes the emission of 
Kelvin waves 



124 



Following this suggestion, Vinen analyzed the energy 
spectrum of the Kelvin-wave cascade by introducing a smoothed length of 
vortex line per unit volume after all the Kelvin waves were removed, and 
considered E-^(k)dk, the energy per unit length of the s moo thed vortex lines 
associated with Kelvin waves in the range k to k + dk 1251 ] . By the dimen- 
sional analysis, E-^ik) was estimated as 

E K (k) = A P K 2 k-\ 



(102) 



with a constant A of order unity. This form is consistent with Eq. ( 11011) . 

Kivotides et al. numerically confirmed the generat ion of Kelvin waves 
through reconnections using the vortex filament model 1261 ] . They also cal- 
culated the energy spectrum E(k) defined by f|44|) for the superfluid velocity, 
and found that E{k) developed approximately a k~ x form. Kivotides sug- 
gested that the fluctuations of the superfluid velocity field were induced by 
the Kelvin waves on the filament, i.e., E(k) ~ E-^(k), and their result was 
consistent with Vinen's analysis of Eq. (11021) . 

The result of Eq. (llOip . Vinen's analysis (I102p . and Kivotides's result 
show the energy spectrum in the quantum regime to be E(k) oc k~ x . How- 
ever, several forms of E(k) have been theoretically predicted by various ap- 
proaches. Kozik and Svistunov analyzed the Kelvin-wave cascade consid- 
ering the three-kelvon scattering process within the weak-turbulence theory 
1291 . Il30j ] and obtain the spectrum E{k) oc k~ 7 / 5 . Nazarenko and Boffetta 
et al. presented a nonlinear differential equation model, pointing out that 
turbulence displays a dual cascade behavior of both the direct energy cascade 
supporting E(k) o c fc~ 7 / 5 and the inverse cascade supporting E(k) oc &T 1 of 



wave action [132l . 1 1 3 II ] . L'vov and Nazarenko also derive the new cascade 
scenario due to the four - Kel vin- wave scattering process showing the energy 
spectrum E{k) oc k~ 5 ^ 3 133]. Yepez et al. suggest the E{k) oc k~ 3 fro m th e 
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numerical result of the GP model as discussed in the previous section 
Boue et al. discuss the discrepancies of these spectrum, especially, the spec- 
tra proposed by Kozik and Svistunov, and L'vov and Nazarenko, and point 
out that it comes from the difference between local (Kozik and Svistunov) 
and nonlocal (L'vov and Nazarenko) theories for Kelvin-wave dynamics jl35 |. 
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They also perform the numerical simulation showing an agreement with the 
nonlocal predictions. Kozik and Svistunov also discuss the importance of 
the symmetry and related Noether's consta nts o f motion in the Kelvin-waves 



when discussing the locality or its absence [145 



3.6.3. Classical-quantum crossover 

As discussed in the previous sections, there are two different types of 
energy spectra in the classical {k < kj) and quantum (fcj < k < regions. 
We have addressed an important question: How do these two energy spectra 
connect to each other at the length scale 17 Although there are several 
theoretical and numerical reports on this region, consistency among these 
works has not yet been obtained. In the analysis of this classical-quantum 
crossover, A = log(//£) appears to be an important parameter. In typical 
4 He experiments, A is about 15. 

Kozik and Svistunov suggested a picture for the crossover region, in which 
the locally induced motion of the vortex lines emerges at the scale of ro ~ 
A 1 / 2 /, and the crossover range is divided into three subranges, r 1 < k < A b , 
K 1 < k < A^, and A" 1 < k < X~\ where A b ~ A 1/4 Z, A c ~ l/K l l\ and 
A* ~ Z/A 1//2 jl37 |. In the first region, r^ 1 < k < A^ 1 , polarized vortex lines 



are organized in bundles and reconnect with other bundles to form Kelvin 
waves with amplitude (k ~ r fc _1 . In the second region, A^ 1 < k < A^ 1 , the 
cascade is supported by nearest neighbor reconnections in a bundle, and (k ~ 
/(A b A;)" 1/2 . In the third range, A c 1 < k < A^ 1 , the cascade is driven by self- 
reconnection of vortex lines, giving Cj, ~ k^ 1 . The Kelvin- wave spectrum 
smoothly connects these ranges. Although they emphasized that the energy 
spectrum E(k) is practically meaningful only in the classical region, we can 
estimate E{k) from their model: E{k) oc k~ 3 in r^ 1 < k < X^, E(k) oc k° 
in A^ 1 < k < X' 1 , and E(k) oc Ar 1 in A,: 1 < k < A^ 1 . 

L'vov et al. suggested another scenario for t he cl assical-quantum crossover: 



bottleneck crossover between the two regions 11361 . For k ~ l' 1 and A > 1, 
the energy of Kelvin waves is much larger than the hydrodynamic energy due 
to superfluid velocity by vortices at the same energy flux. As a result, there 
is a bottleneck energy accumulation around k ~ l" 1 and the energy spec- 
trum becomes E(k) oc k 2 for superfluid velocity and E^(k) oc k° for Kelvin 
waves, followed by equipartition of the hydrodynamic energy and the energy 
of Kelvin waves, respectively. This scenario is completely different from that 
suggested by Kozik and Svistunov because there is no energy stagnation in 
the model by Kozik and Svistunov. 
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Figure 24: Incompressible energy spectra plotted vs. kl. The simulation results (the 
same symbols as in Fig. |22|) and the model by L'vov et al. for A = 10, 30, 100 (dashed 
curves) are brought together to the theoretical (solid) curve with A = 2 by superposing 
the Kolmogorov spectrum (for both simulations and model) and plateau regions (only 
for simulations. The dot-dashed lines show different scaling asymptotes. [Sasa, Kano, 
Machida, L'vov, Rudenko, and Tsubota: Phys. Rev. B 84 (2011) 054525, reproduced 
with permission. Copyright 2011 the American Physical Society] 



As discussed in Sec. 13.6.11 the large-scale numerical simulation of the 
GP model quantitatively supports the existence of a bottleneck effect in the 



crossover region 142| . Figure [2H shows a comparison between the numerical 
result and the theoretical prediction for different A. (For the sake of a better 
comparison we replotted the simulation data, bringing them all together to 
the model by L'vov et al. with A = 2 by superposing the Kolmogorov 
spectrum and plateau regions.) However, A is too small to predict the exact 
mechanism and the comparison shown in Fig. 1241 remains problematic. More 
theoretical studies and numerical and laboratory experiments are required to 
fully understand the vortex dynamics in the scale crossover region. 

4. Quantum hydrodynamic instability in two-component Bose— Einstein 
condensates 

Hydrod y na mic i nstability is of fundamental importance in classical fluid 



dynamics |146l . 11471 ] . The instability causes characteristic wavy patterns de- 
veloping into turbulent flows from a basic flow after complex dynamics of 
eddies. Because of the universal applicability of the theory, hydrodynamic 
instability appears in different kinds of fluids. SuperfTuids are no exception. 
Although superfiuid dynamics are described by hydrodynamic equations 
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similar to those in classical fluid dynamics, the macroscopic quantum effects, 
i.e., superfluidity and vortex quantization, can cause significant differences 
between the hydrodynamic instabilities in quantum and classical fluid sys- 
tems. Superfluidity, i.e. flow without friction, enables us to study different 
hydrodynamic instability arising from frictionless basic flows. Such an in- 
stability has no classical counterparts because frictionless flows cannot be 
achieved in classical fluid systems because of the presence of viscosity. On 
the other hand, the appearance of a quantized vortex should cause a sig- 
nificant difference between the quantum and classical fluids at least in the 
nonlinear stage of instability where vortices are relevant. Therefore, quan- 
tum effects cause differences in both linear stability and nonlinear dynamics 
between quantum and classical hydrodynamic instability. 

In this section, we discuss hydrodynamic instability in quantum flu- 
ids, namely, quantum hydrodynamic instability, especially in two-component 
BECs. Historically, the study of quantum hydrodynamic instability has been 
developed in helium superfluid systems. First, we shall briefly review hy- 
drodynamic instability in superfluid systems focusing on helium superfluid 
systems. Then, we introduce hydrodynamic instability in two-component 
BECs. 



4-1. Hydrodynamic instabilities in superfluid systems 

The Landau instability 148] is the fundamental mechanism for the linear 
stability of frictionless flows. The Landau instability occurs when frictionless 
states become unstable when the superflow velocity relative to the external 
environment, such as the container wall or the thermal excitations dragged 
by the wall, exceeds a critical value. This instability is a thermodynamic 
instability, where elementary excitations with negative energy are sponta- 
neously amplified in the relaxational process, decreasing the thermodynamic 
energy of the system. 

The most important example of quantum hydrodynamic instability is the 
instability of thermal counterflow in superfluid 4 He. The thermal counterflow 



instability has been studied in parallel with QT in superfluid 4 He [22j, [23, 



24j . |25| | as was introduced in Sec. 13.31 The thermal counterflow instability 
occurs when the relative velocity between the normal fluid and superfluid 
components exceeds a critical value. Then, remnant vortices attached to the 
container wall are stretched by the mutual friction, and the stretched vortices 
repeat reconnections, leading to QT. 
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The stretch of remnant vortices in the thermal counterflow instability may 
be understood more fundamentally as an instability of single quantized vor- 
tices, the Kelvin- wave in stability, sometimes called the Donnelly-Glaberson 
instability 1491 . Il50[ Il51| . Kelvin- wave instability can occur when there is a 



relative helical flow along a vortex line between normal fluid and superfluid 
components, which is realized by injecting a heat current along the quantized 
vortices. The instability leads to amplification of the Kelvin waves and helical 
deformations of the vortex line. When the Kelvin-wave instability is induced 
in r otating su perfluids with vortex lattices, the instability can develop into 



QT [152l . Il53| . The Kelvin- wave instability and its counterpart have been 
discussed for other superfluids and superconductors. The instability can oc- 
cur in superfluid 3 He-B in a rotating cylinder, where the rotating counterflow 
of a vortex -free superfluid component and a rotating normal fluid component 
is realized 
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In atomic BECs, it was proposed that the Kelvin-wave in- 
stability can be induced as a spontaneous exci tatio n of a kelvon, a quantum 
of a Kelvin wave due to the Landau instability 155] . The Kelvin- wave insta- 



bility has also been discussed for quantized vortices in a rotating neutron star 
156l . Il57| . In type-II superconductors, the counte rpart of the Kelvin- wave 



instability is the spiral- vortex expansion instability 158| , where a flux vortex 
becomes unstable against the growth of helical perturbations in the presence 
of a sufficiently large current density applied parallel to its axis. 

The above examples of quantum hydrodynamic instability are phenom- 
ena that do not occur in classical fluids. Also of interest is the study of 
the quantum counterpart of classical hydrodynamic instability. The Kelvin- 
Helmholtz instability, one of the most fundamental instabilities in classical 
fluid dynamics, was first studied experimentally in a superflu id sy stem at 



the interface between the A and B phases of superfluid He [1591 ] . When 
the relative velocity between the two phases exceeds a critical value, the 
penetration of quantized vortices across the interface was detected by count- 
ing the number of vortices by NMR before and after the instability. The 
Kelvin-Helmholtz instability in superfluid systems has been dis cussed on 
different kinds of interfaces, such as n orm al-superfluid interfaces [lloL Il60| . 



nuclear-nuclear superfluid interfaces [1611 ]. and superfluid-superfluid inter- 



lr The instability observed in the experiment [159j | is the thermodynamic instability 
triggered by the Landau instability of ripplons of the interface between the A and B phases. 
Strictly speaking, this instability is not the counterpart of the classical Kelvin-Helmholtz 
instability, which is the dynamic instability of ripplons. See the following section. 
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face in atomic two-component BECs 162 



Recently, there has been growing interest in hydrodynamic instability in 
atomic BECs. The unique dynamics due to quantized vortices have been re- 
ported on the quantum counterp arts of classical hydrodynamic in stabilities , 
the Kelvin-Helmhol tz in stability 162| , Rayleigh-Taylor in stab ility 1631 . Il64 



Strouhal instability [165], Richtmyer-Meshkov instability 166], and Plateau 



Rayleigh (capillary) instability |l67j , among others. There are several merits 
in considering an atomic BEC system to study quantum hydrodynamic insta- 
bility. The most important advantage of such a system over other superfluid 
systems is that we can visualize directly the whole time development of the 
order parameters in the instability dynamics, from the linear stage to the 



168 . In addition, since 



nonlinear development including vortex nucleation 
the system is less dissipative at ultra low temperatures, we can observe the 
hydrodynamic instabilities, which are obscured by dominant thermodynamic 
instability in dissipative systems. The theoretical advantage of this system 
is that we can predict quantitatively the detailed dynamics of quantum hy- 
drodynamic instability within the mean field approximation. 

Multi-component atomic BECs provides an ideal ground to study novel 
hydrodynamic phenomena of multi-superfluid systems. Two-com pone nt atomic 
BECs are the simplest systems of multi-component superfluids (169| . which 
can be created in cold-atom systems with multiple hyperfine spin states or 
a mixture of different atomic species. Recent experimental advances enable 
us to study a variety of superfluid dynamics in two-component BECs in a 
more controllable manner. The intra- and inter-compon ent interaction s can 
be tuned with the help of the the Feshbach resonance [173, [ml, Eel; and 
the external potentials are controllable independently on both components 
by utilizing the difference between the Zeeman shifts of the two components. 
Recently, hydrodynamic instability of counterflows in miscible two superflu- 
ids, called countersuperflow instability, was observed for the first time in 



two-component BECs by Hamner et. al. [1731 ]. Coutersuperflow instability 
is unique to mult i- component superfluid systems. It was also suggested that 
countersupe rflow instability can develop into a binary QT composed of two 
superfluids jl74| . 

In the following subsection, we introduce our recent work on quantum hy- 
drodynamic instability in two-component BECs. The next subsection is de- 
voted to introducing the hydrodynamic formalism for two-component BECs 
to make it easier to understand the subsequent subsections. Then w e develop 



the discussion into concrete problems, countersuperflow instability [I74j . 1175 
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and quantum Kelvin- Helmholtz instability [162l . Il76j . which are the most 
fundamental hydrodynamic instabilities in two-component BECs. 



4-2. Linear stability and hydrodynamic formalism 

Let us start with the Lagrangian of two-component BECs 



h 
i— 

2 



JdVj2 ~ VAVj) - K, (103) 



where the index j refers to the jth component. If we take into account the 
motion of an external environment with the translational velocity V ex and 
the ang ular velocity fl ex , the thermodynamic energy K is generally written 



as 



177 



K 



J dV[K- (V ex + Q ex x r) ■ J] , (104) 



where K. is the thermodynamic energy density for V ex = Q ex = 0, and 
J = Yi^j^j^^i ~ ^i^^j) is tne tota l momentum density of the two 
components. The energy density K is written as K = fC\ + IC2 with 



K. 



h 2 



3 



I v^f + (Uj - ^Ol^l 2 + 9jk\*3\ 2 \*k\ 2 , (105) 



2m,- ' •" v J J " •" 2 
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where mj, Uj(r), and fij are the particle mass, the external potential, and 
the chemical potential of the jth component, respectively. The inter- and 
intracomponent interaction constants g^ have the form g^ = 2irh 2 ajk(mJ 1 + 
m^ 1 ), where is the s-wave scattering length between the jth and fcth 
components. For simplicity, the environment is supposed to be at rest in the 
laboratory frame throughout our discussion, namely V ex = Q ex = 0. 
From the Lagrangian (I103p . one obtain the coupled GP equations, 



ihd t fyj 



h 2 



2m,- 

3 k 



(106) 



Steady superflows are realized as a stationary solution *&j(r,t) = &j(r) of 
the GP equations (j!06p . For example, without external potential Uj = 0, 
the GP equations (11061) have the stationary solutions <&j(r) oc ^ im i v i- r l h of 
steady uniform superflows with arbitrary velocity V ~. 
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To clarify the relation to fluid dynamics, we introduce a hydrodynamic 
formalism. By inserting ^ = fje %Bj into the Lagrangian density C, we obtain 

( v v^) 2 + ( u i - ^j) n i + \ ^dikUkUj L (107) 

k J 



h 2 



where rij = f? and Vj = ^V0j is the density and the superfluid velocity 
of the jth component, respectively The variation with respect to fj and 9j 
yields a set of hydrodynamic equations 



dtnj + V • (rijVj) = 0, (10 
rrijdtVj = — V 



^v] ■< h -r J - ,!>;_. (109) 



where q,j = _ ^:(V 2 /j) / /j is the quantum pressure term, and fx^ = J2 k g 3 k n k 
is the hydrostatic chemical potential, which is named after the hydrostatic 
pressure in fluid dynamics. These hydrodynamic equations are analogs of 
those of multi-phase fluids in classical fluid dynamics. The first equation 
represents the conservation law of density rij and the second has a similar 
form to the Euler equation of irrotational flows. If we neglect the inter- 
component interaction, g 12 = 0, the term V/ij reduces to ^-Vp^ with the 
hydrostatic pressure 

V) Uu"] (HO) 

of the jth component. In a stationary state, the second equation reduces to 

-^-v) + qj + Uj + n) = Hj = const., (Ill) 

which is the counterpart of the Bernoulli theorem. 

The interaction between different components comes from the term g^k (j 7^ 
k) in the hydrostatic chemical potential \i h j. The force density fj k on the jth 
component by the fcth component is written as 

fjk = -9jknjVn k . (112) 
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The inter-component force Fj k = f dV f ]k obeys the principle of action 
and reaction, Fj k = —F k j. The force vanishes for uniform density profiles 
rij = const., which makes it possible to realize a stationary relative flow 
v\ 7^ i>2 between different components, namely, a countersuperflow state. 

Linear stability analysis can be done in a manner similar to that in clas- 
sical fluid dynamics. To perform the linear stability analysis, we introduce 
small perturbations 

rij(r,t) = fij(r) + 8nj(r,t), (113) 
O s (r,t) = ^(r) + ^(r,t), (114) 
Vj {r,t) = v 3 (r)+5v 3 (r,t), (115) 

where Vj = ^W,- and Svj = £rV68j. By linearizing Eqs. f TT08|) and ( TT09]) . 
we obtain the linearized equations 

d t 5rij + V • (nj5vj + 5njVj) = (H6) 
rrijd t 8vj = — V [rrijVj ■ 5vj + 5qj + 8jJ^\ , (H7) 

where 

h 2 

5 Qj = - — —[^.(V 2 ^) + n^Vnj) ■ V - (Vfij) 2 - n^V^rij (118) 

and 5^ = Y2k 9j^ n k- These equations determine the linear stability of the 
stationary state. 

As the first example we consider the linear stability of uniform superflows 
without relative velocity: v = V\ = v 2 = const., Vj(r) = const, and fij(r) = 
const. The perturbations may be written as 5rij oc cos(q • r — ut), 59 j oc 
sin(qr • r — ut) with the wave number q and the frequency u. The linearized 
equations (11161) and (I117P are reduced to the eigenvalue equations for the 
eigenvalue u = u — v ■ q, and we obtain the dispersion 

co = v ■ q ± \J\{^1 + ul) ± l -^o\ - u 2 ) 2 + 4cf 2 q\ (119) 

with to 2 = c 2 ,q 2 + -r-rQ 4 and c 2 h = \/q% " J,lt . The first term comes from the 

Doppler shift due to the background superflow v and the second term refers 
to the eigenvalue uq. 
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The thermodynamic stability of superflows or the Landau instability is 
investigated from the deviation SK of the thermodynamic energy K due to 
perturbations. By using the equations f II 1 6 j) and (11171) . the deviation energy 
SK is written as 

SK = - / dV^2 (SOjdtSrij - dnfitdOj) . (120) 
2 J . 

The deviation energy SK due to the perturbation is proportional to the 
frequency u (II 19ft . Since u changes its sign by increasing the velocity v — \v\ 
with the eigenvalue uq fixed, SK becomes negative for a perturbation when 
v exceeds a critical velocity vl, called the Landau critical velocity. Then 
the perturbation is amplified to decrease the thermodynamic energy of the 
system. From the dispersion (I119p . the Landau critical velocity is given by 




The frequency can take a complex value, where the perturbations are ex- 
ponentially amplified as oc e at with a = Im u > 0. Then the system is called 
dynamically unstable. The dynamic instability occurs when g\ 2 > g\\g22 for 
the dispersion (11190 . The instability means that homogeneous condensates 
are unstable; if g±2 > -y/gnfi^, the two components will undergo phase separa- 
tions due to the strong inter-component interaction, while if —gu > ^gng22, 
the condensates are unstable to the formation of a denser droplet containing 
both components. These dynamic instabilities are classified as hydrostatic 
instability rather than hydrodynamic instability in the sense that they occur 
without superflow, Uj = 0. 

Dynamic instability is purely an internal instability in isolated systems 
without external environments, where the particle numbers Nj = J dVrij 
and the energy E = K + ^ ■ (J,jNj are conserved. On the other hand, the 
Landau instability in this case is a hydrodynamic instability induced by the 
frictional dissipation between the condensates and the external environment. 
If the system under consideration is dissipative, rather than dynamically 
unstable, the Landau instability is dominant on hydrodynamic instability. 
The frictional dissipation can be small and thus dynamic instability can be 
dominant in atomic BECs when condensates are trapped in 'a smoothed-wall 
container' made by electromagnetic fields, at low temperatures, where there 
is a small amount of the normal fluid component. 
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In the following subsections, we discuss the countersuperflow instability 
and the Kelvin-Helmholtz instability, which are fundamental instabilities 
in the presence of relative velocity between different components. These 
instabilities belong to the dynamic instability attributed to the relative flow 
of two condensates, although the Landau instability occurs due to the relative 
motion between the environment and condensates. 



4-3. Countersuperflow instability in miscible two -component BECs 

It has already been mentioned that countersuperflow states can be re- 
alized as a basic flow in two-component BECs. It is interesting to compare 
countersuperflow with thermal counterflow of superfluid 4 He. Both states are 
a counterflow of two miscible fluid components, which has no analog in classi- 
cal fluid dynamics. The latter is a counterflow of normal fluid and superfluid 
components, which becomes thermodynamically unstable when the relative 
velocity exceeds a critical value. The former is a counterflow of two super- 
fluids. At first sight, it seems that the countersuperflow state can be stable 
for arbitrary relative velocities since each component is a superfluid by itself. 
However, the countersuperflows in two-comp onen t BE Cs become dynami- 



cally unstable over a critical relative velocity [1781 . Il79l | . In this subsection, 
we discuss the linear stability and nonlinear dynamics of countersuperflow 
instability in two-component BECs. 

4-3.1. Linear stability of countersuperflows 

Let us consider uniform counterflows in two-component BECs in an iso- 
lated homogeneous system, rij = fij = const, and vj = Vj = const, with 
relative velocity vr = v 2 — V\ 7^ 0. Because of the Galilean invariance, 
we can neglect the translational motion of the whole system without loss 
of generality, (Pi + P 2 ) = M\V\ + M 2 v 2 = with the total momentum 
Pj = m j J dVnjVj and the total mass Mj = mjNj of the jth component. 

The linear stability of the countersuperflows is evaluated using the cou- 
pled linearized equations (j!16p and All Tj) . The problem is reduced to solving 
the equations 

(d t + v* ■ VfSrii = -^V 2 (5/i 7 - + SqA. (122) 

rrij 

By substituting 5rij oc cos(q • r — ut) into Eq. (I122p . we obtain 

[(u - v ± ■ q) 2 - ul][(u - v 2 ■ q) 2 - u 2 ] = c\ 2 q 4 . (123) 
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Although the stability is investigated by solving this quartic equation, the 
eigenvalue has a complicated form in general. However, for a symmetric case 
mn = m 22 = m, rii = n,2 = n, and gw = #22 = g, which is a reason able 



approximation, e.g., for two-component BECs of 87 Rb atoms 1721 . |l73|, the 
dispersion relation is reduced to a simple form 



s' 2 = \q' A + q' 2 + ± sj C-q* + qp)tfV$ + (124) 

where e' = frw/gn, q' = q£ with £ = h/ ^/mgn, V R = \V' R \ = \v R \/c with 
c = a/ gn/m, and 7 = gn/g- Here q' 2 = gf 2 + with q'» = \q' ■ V' R /V R \ 
and q' ± > 0. By comparing the last term with the sum of the other terms on 
the right hand side in Eq. (I124p . the condition Im e' ^ for the dynamic 
instability is found to be 



+ q ,2 {l _ 7 ) < < ^/l q /4 + q ,2 {1 + 7 ). (125) 

It is reasonable to expect that this inequality is never satisfied without the 
inter-component interaction (7 = 0). 

Figure [251 shows the phase diagram of the countersuperflow instability for 
7 = 0.9. The unstable region is characterized by the lower and upper critical 
velocities 

V' ± = 2^(1 ±|7l)- (126) 

The unstable region appears when V' R exceeds the lower critical velocity V'_- 
The distribution of the unstable modes, which have larger values of |Im e'\, 
depends on the relative velocity. For sufficiently large relative velocity with 
V R > V+, the cross section of the unstable region has a crescent-like form 
[Figs. [25] (a)]. In this case, we find that the unstable mode with large values 
of |Im e'\ are distributed in the region with higher wave number q' ± . On 
the other hand, if V' R decreases below V+, the unstable region is broadly 
distributed around q' ± = [Fig. [251 (b)]. 

4-3.2. Nonlinear development of countersuperflow instability 

The distribution of the unstable modes in the wave number space strongly 
affects the nonlinear dynamics of the vortex nucleation after the linear am- 
plification of the modes. Figure [26] shows a typical time development of 
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Figure 25: Phase diagram of the countersuperfiow instability for 7 = 0.9. The dashed 
lines represent (a) V' R = 4.71 and (b) 2.36 and the solid lines show the critical velocities 
V' ± = 2^/1 ± 7. The two right-hand plots show the cross-section surfaces of V' R — 4.71 
and 2.36 of the phase diagrams. [Ishino, Tsubota and Takeuchi: Phys. Rev. A 83 (2011) 
063602, reproduced with permission. Copyright 2011 the American Physical Society] 

the countersuperfiow instability, obtained by numerical simulation of the GP 
equations (11061) . The numerical simulations were done in a three-dimensional 
box under periodic boundary conditions with the parameter settings ni\ = 
m 2 = m, g u = g 22 = g, n x = n 2 = n, v x = -v 2 , = mv 2 R /8 + (g + g 12 )n, 
7 = 0.9, and V R = 4.71. We add a small amount of white noise in the initial 
countersuperfiow state to trigger the instability Because of the symmetric 
parameters between the two components, the nonlinear dynamics are similar 
for both components. 

After the exponential amplification of the unstable modes, disk-shaped 
low-density regions appear in both components, which face in the direction 
parallel to the initial relative velocity [Fig. [261( b)] . and then vortex rings are 
nucleated inside the regions [Fig. l26T c)]. Since the size of the vortex rings 
is similar to those of the low density regions, the vortex distribution is char- 
acterized by the density pattern emerging after the onset of the instability. 
Thus, if the instability is so strong that the density pattern grows soon into 
vortex rings, the vortex number density immediately after the instability is 
estimated by the wavelength of the most unstable modes in the phase dia- 
gram. For the limit of large relative velocity V' R ^> VI , the unstable region 
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Component 1 (vortex core and density isusurface) 








Component 2 (vortex core and density isosurface) 
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(a) f = (b)i'=12.2 (c)t'=12.S (d)t'=13.3 (e)f'=I3.8 (f) f = 20.0 



Figure 26: Nonlinear dynamics of vortex cores in countersuperfiow instability. The 
top and middle panels show the vortex dynamics in the first and second components, 
respectively. The surface plots represent the low density isosurface of |Wi| 2 /n = 0.1 (top) 
and |* 2 | 2 /ri = 0.1 (middle). [Ishino, Tsubota and Takeuchi: Phys. Rev. A 83 (2011) 
063602, reproduced with permission. Copyright 2011 the American Physical Society] 



determined by the inequality ( 1125ft is reduced to (ql — \V^) 2 + q'± = \V R 2 , 
and then the vortex line density is estimated to be ~ v 2 r /k 2 . 

For the case of small relative velocity, e.g. V'_ < V' R < V+, the size of the 
nucleated vortex rings becomes large because the unstable modes with larger 
imaginary part are distributed mainly around q± = 0. If the perpendicular 
wavelength 2ir/q±_ of the unstable modes is similar to or larger than the sys- 
tem size normal to the relative flow, the instability does not cause nucleation 
of vortex rings, but rather makes vortex lines acros s the system or distorted 



stripe patterns, as observed by Hamner et al. |173 . 

Typical dynamics after the vortex ring nucleation are demonstrated in 
Figs. l26T d)-(f). The vortex rings propagate along the initial relative ve- 
locity in the opposite direction between the two components, and the ring 
size increases with time. When vortex rings come close to each other, the 
rings are distorted by the interaction and make reconnections [Fig. [261(e)]. 
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The distortion and the reconnection mainly occur between the vortices in the 
same components because the vortex-vortex interaction be tween the same 
component is larger than that between different components 180| ; reconnec- 
tion does not occur between vortices in the different components. Quantized 
vortices in both components become tangled due to the distortion, forming 
binary QT [Fig. E2T)]. 




Figure 27: Time evolution of the total vortex line density I = l'£~ 2 and the velocity 
V x .j — v' X j^/T of the jth component along the initial relative velocity. The time and 
length are scaled by £ = h/ ^/mgn and r = %/ fj,. 

The time development of the countersuperflow instability can be consid- 
ered as the frictional relaxation of the relative motion of the two condensates. 
In general, a frictional force between two interacting objects causes a decrease 
in their relative motion and their kinetic energy is dissipated in various forms 
of energy, i.e., the internal energy such as heat. Here, we define the macro- 
scopic kinetic energy Ekin of the two condensates as 



E, 



kin 



p2 
r 2 



+ 



p2 
r 2 



2Mi 2M 



2 



R 



(127) 

Pi /Mi and the 



with the macroscopic relative velocity V r = P2/M2 
macroscopic reduced mass M' = (Mf + M% The corresponding 'in- 

ternal energy' E int is defined by subtracting E kin from the total energy K\ 
Ei n t = K — Ekin- Then, the frictional force Fr between the two compo- 
nents can be defined from the time variation of the relative kinetic energy 
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E R = \M'V\. 




F R -V 



(128) 



with the frictional force 



Fr EE M'-Vr. 



(129) 



Since M' is the conserved quantity, Fr is proportional to the time variation 
of the relative velocity V r. 

Figure I2T1 shows the time evolutions of the averaged velocity Vj = Pj/Mj 
of the jth component and the sum of the vortex line density of the two compo- 
nents. Throughout the development, the velocity Vj is almost perpendicular 
to the initial velocity Vj \\ x with the unit vector x along the x-axis, and only 
the parallel component V x j = Vj ■ x is plotted in Fig. [271 From the relation 
(11291) . the reduction rate of \V X> 2 — V Xj \\ is proportional to the frictional force 
between the two components. 

The friction between the two components is small in the linear stage of 
the instability, but grows drastically as the vortex line density increases. The 
decrease in the kinetic energy Er by the friction is offset with the increase in 
the 'internal energy' by nucleating and expanding vortex rings. On the other 
hand, the vortex reconnection suppresses the friction since the distortion 
of the vortex ring configuration due to the reconnection disturbs the free 
expansion of the vortex rings. In addition, the vortex line length can be 
decreased after vortex reconnections since some of the energy is dissipated 
for phonon emission. The vortex line density increases to a maximum value 
when the two effects, the vortex ring expansion and the vortex reconnection, 
are balanced in the vortex-tangled state. The frictional relaxation continues 
but its rate decreases in the tangled state. The total length starts to decrease 
when the relative velocity becomes almost zero, and then the binary QT will 
decay. 

The dynamics of the turbulence transition in countersuperflow instability 
is similar to that in thermal counterflow instability. Recall that in the ther- 
mal counterflow instability, quantized vortices are stretched by the mutual 
friction between the superffuid and normal fluid components. The steady 
QT developed from the thermal counterflow instability is anisotropic since 
the relative velocity between the two components is sustained externally by 
applying a temperature gradient through the system. On the other hand, 
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the counterpart of the mutual friction is caused by the friction between the 
two condensates in the countersuperflow system. Quasi-steady turbulence 
is realized temporarily when the vortex line density approaches the maxi- 
mum value. It is numerically shown that the maximum vortex line density is 
proportional to the square of the initial relative velocity, similarly to the rela- 
tion of Eq. (T2Tj) . but the vortex tangle can be isotropic w hen t he momentum 



exchange is completed and the relative velocity vanishes [175 



4-3.3. Conclusion 

Countersuperflow states in miscible two-component BECs become dy- 
namically unstable when the relative velocity between the different compo- 
nents exceeds a critical value. The countersuperflow instability causes vortex 
nucleation and stretching vortices leading to isotropic binary QT. The time 
development of the countersuperflow instability is interpreted as frictional 
reduction of the relative motion of two condensates. These phenomena are 
interesting in two senses. One is that the relative superflows decay due to 
the mutual friction between two superfluids, each of which consists of the 
'frictionless' superfluid component by itself. The other is that the QT of 
mult i- component BECs can be realized from the countersuperflow instabil- 
ity. These phenomena can be observed with current experimental techniques 



1731 . |181| if the size of the condensates is sufficiently large in the direction 
perpendicular to the relative velocity. We hope that these phenomena will 
be observed in future experiments. 

4-4- Kelvin-Helmholtz instability in immiscible two- component BECs 

In classical fluids, the Kelvin-Helmholtz instability can occur when there 
is a sufficient velocity di fference a cross the interface between two fluids with 



different mass densities [1461 . 11471 ] . A vortex sheet exists along the interface 
due to the velocity difference and the instability induces exponential amplifi- 
cation of the oscillating modes of the vortex sheet. The instability typically 
develops into roll-up patterns of the interface in the nonlinear stage. The 
Kelvin-Helmholtz instability is one of the most fundamental hydrodynamic 
instabilities in fluid dynamics, related to several familiar phenomena such as 
wind-generated ocean waves, flapping flags, billow clouds, and sand dunes. 
In this subsection, we discuss hydrodynamic instability in phase-separated 
condensates, where the inter-component interaction parameter satisfies the 
immiscible condition g±2 > y/gwg-n- We shall show that, in the presence of 
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the relative velocity between the phase-separated condensates, the instability 
is related to the Kelvin-Helmholtz instability in classical fluid dynamics. 

4-4- 1- Linear stability of a flat interface 





iii 
Interface layer 
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z/Z 

Figure 28: Profiles of the order parameter amplitudes fj = ^Jn] of phase separated two- 
component BECs under external potentials Uj(z) — GjZ for the jth component. The real 
part of the excitation function Srpj = Uj(z) — Vj{z)* of a typical interface mode is plotted 
with broken lines. The parameters are set as /i = ju,j — -?r^, m — rrij, g — gjj, gn = 10g, 
and Gj — — (— l) J G/2 with G — 0.02/i/^. The amplitude and the length are scaled by fi/g 
and ^ = H/y/mfi, respectively. [Takeuchi, Suzuki, Kasamatsu, Saito and Tsubota: Phys. 
Rev. B 81 (2010) 094517, reproduced with permission. Copyright 2010 the American 
Physical Society.] 

We consider a flat interface in phase-separated condensates. The interface 
is located at z — 0, where fx > f 2 for z < and / 2 < fi for z > as shown 
in Fig. |28j The position of the interface is stabilized at z = under the 
small potential gradient d z Uj = Gj = const, with G\ > and G 2 < 0. The 
stationary solution has a form <&j(r) = fj(z)e m i Vrr ' h , where the superfluid 
velocity Vj is parallel to the interface. The interface layer can be defined 
as the region sandwiched between the regions m ~ nj and n 2 ~ n 2 , where 
nj(z) = (fij — Uj (z) — I gjj is the bulk density obtained by neglecting the 
quantum pressure and the density of the different components = (k ^ j) 
in the Bernoulli equations f illip . The thickness of the layer decreases when 
the inter-component interaction becomes large. For simplicity, we consider 
strong phase separation with sufficiently large gi 2 , where the thickness is 
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minimized to the order of ~ fx/ ymjgjjnj (0), the healing length of a 
single condensate. 

The linear stability of the stationary states $j is investigated by lin- 
earizing the GP equations with respect to a collective excitation S^j(r,t) = 
*&j(r, t) — $j(r). An oscillating perturbation of frequency u is conventionally 
described with the Bogoliubov formalism 



V " n { Uj (z 



iqr—iuit 



[Vj(z}e 



iq-r—iu)t~\ * 



} 



(130) 



where the wave number q is parallel to the interface. The functions Uj and 
Vj obey the reduced Bogoliubov-de Gennes (BdG) equations, 



ti 2 ( m 7 - _ 



h 2 d 2 

2rrij dz 2 



2m 



+ ^2djk (fkUj + fjfkUk ~ fjfkVk) = fctt 
k 

h 2 d 2 



j- 



(131) 



m 



-2.V- 

h 3 



2m j dz 2 



^2 g i k (fk v j + fjf^k ~ fjfkUk) = hwvj- (132) 



These equations determine the linear stability of the stationary states. 

Excitations in the phase-separated states are generally classified into two 
types. One is bulk modes such as phonons, which can propagate in the bulk 
far from the interface. The other is localized modes, which disturb the order 
parameters only locally around the interface and decay exponentially in the 
bulk (see Fig. 128]) . The simplest example of localized modes is a transverse 
shift of the interface in the z direction. In a manner similar to the linear 
stability analysis of the Kelvin-Helmholtz instability in hydrodynamics, we 
consider only oscillations of such interface modes, neglecting the bulk modes 
and the internal structure of the interface. 

If the thickness of the interface is neglected and the interface position is 
represented by the single- valued function z = rj(x,y,t), the interface modes 
may be approximately described by the effective Lagrangian 



^eff 



dxdy 



dzV\ + / dzVi — Q.S 



(133) 
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where we used 



v, = -hn 3 d t e 3 - '-^vj + ^/iV a /i - (Uj - N )m - \g S fi#M) 

S = y/l + (d x ri) 2 + {d y ri) 2 (135) 

and the interface tension coefficient a was introduced. Let us consider a flat 
interface r\ = in a stationary state in a homogeneous system along the x 
and y axes. For a small perturbation, we obtain the equation of motion for 
the interface position r], 

Vi(v) - V 2 {r]) + aV 2 r] = 0. (136) 

This equation is an analogue of the Bernoulli theorem on the interface. In 
the stationary state rj — 0, the term V\ is reduced to the hydrostatic pressure 
Pj f lllUp of the jth component at the interface; we then obtain p\ = p\. 

We can employ the kinematic boundary condition on the interface, similar 
to the discussion for hydrodynamics, 

(Sv^^dtrj + Vj-Vr). (137) 

Based on the assumption that the density perturbation 5rij at z ~ is caused 
by the local transverse shift of the density profile, we can write the density 
perturbation along the interface as 

5rij = —r]d z rij. (138) 

To obtain the dispersion of the interface modes, we assume the localized 
perturbations with a form 77 oc sm(q-r—cut), 8n 3 - oc e _< - _1 ^ aj2: sin(q-r— cot) and 
58 j oc e _< - _1 ^ ajZ cos(q • r — cot). We obtain a 3 = q = \q\ in the approximation 
neglecting the quantum pressure term, d z n 3 m d z nj = —(—iyG/2gjj, by 
linearizing Eqs. (jll3j) . f )137p . and (I138p with respect to the perturbations. It 
is straightforward to derive the dispersion relation for the interface mode, 

u(q) = q-v G ± lq(F + aq*) - -^^(v R ■ qf, (139) 

VPl + Pi V Pl + P2 

where pj = rrijnj(0), F = Ginj(0) — G , 2^2 1 (0), Vr = v 2 — t>i, and v G = 
P1V1+P2V2 

P1+P2 
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Dynamic instability occurs when the imaginary part Im u becomes nonzero 
for F+ ^ q < -£±£L- V 2 R with vr = \vr\. This instability is the counterpart of the 
Kelvin-Helmholtz instability in classical fluid dynamics. The critical relative 
velocity Vd for the dynamic instability is given by 



P1P2 

Note that the interface is hydrostatically unstable without the relative ve- 
locity for F < 0. This is an analog of the Rayleigh-Taylor instability, where 
a heavier fluid is located above a lighter fluid and the interface between the 
two fluids is unbalanced under gravity. 

The Landau instability is evaluated from Eq. (I139p . We obtain the 
Landau critical velocity Vl for the velocity vq = \Vg\ from the condition 
uj = 0, 



Vl = -==j2VF^--^-vf l . (140) 

VPl + P2 V Pi + Pi 

The Landau critical velocity Vl depends on the relative velocity vr. When 
vr > Vd, the system becomes thermodynamically unstable for an arbitrarily 
small velocity vq- On the other hand, even when vq > Vl, the system can 
be still dynamically stable with Im u — 0. Therefore, the Landau instability 
can occur in general before the onset of dynamic instability in a dissipative 
system. The Landau instability has been e xperiment ally observed at the 



interface between the A and B phases of 3 He 1591 . Il82 . 

Figure [221 shows the phase diagram of the dynamic instability (DI) and 
the Landau instability (LI), obtained by the dispersion (I139p . Here, we 
considered q \\ Vr and V\ — with the parameters described in the caption 
of Fig. |28l The results are compared with those obtained by the direct 
numerical computations of the BdG equations (11311) and (11321) . The analytic 
results are in good agreement with the numerical results. Note that the 
approximation demonstrated here is not applied to the perturbation with 
q > l/£ since the interface thickness ~ £ is neglected in this model. The 
difference for small q comes from the inadequate treatment of the density and 
phase perturbation when the penetration depth a -1 = q^ 1 of the interface 
modes becomes comparable to the system size. 

4-4-2. Nonlinear development of dynamic Kelvin-Helmholtz instability 

The nonlinear time development of the dynamic and Landau instabilities 
is investigated numerically. We shall demonstrate typical developments of 



84 




0.4 

0.4 0.8 

Figure 29: Phase diagram of the dynamic instability (DI) and the Landau instability 
(DI) for V\ = and V% < 0. The curves are the boundaries of the DI and LI regions 
obtained from numerical calculation of the BdG equations (solid curves) and from the 
dispersion of Eq. (I130[) (broken curves). The parameters are the same as th ose i n Fig. 
|2"51 The interface tension is calculated as a = 0.886/i 2 £/ g according to Ref. [l83| . The 
relative velocity vr and the wave number q are scaled by c = \J gp/m 2 and £ = h/ ' ^fgp 
with p — pj — m/i/g. [Takeuchi, Suzuki, Kasamatsu, Saito and Tsubota: Phys. Rev. B 
81 (2010) 094517, reproduced with permission. Copyright 2010 the American Physical 
Society] 
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instabilities in quasi-two-dimensional systems neglecting the y coordinate. 

We first discuss the nonlinear time development of the dynamic insta- 
bility (dynamic Kelvin-Helmholtz instability) for vr > Vp. The nonlinear 
dynamics of the Kelvin-Helmholtz instability are obtained by numerically 
solving the GP equations f ll06[) with a quasi-two-dimensional system peri- 
odic along the relative velocity. Figure EDI shows a typical time development 
of the density difference n\ — n 2 in the dynamic Kelvin-Helmholtz instability. 
The vorticity u v and the mass current velocity v, 

co v = rot V, V = Jl+j2 (141) 

m\n\ + m 2 n 2 

are useful for understanding the phenomena. When there is a velocity differ- 
ence between components across the interface, the vorticity co v is distributed 
along the interface. A quantized vortex in the bulk far from the interface has 
a singular peak at its core in the vorticity distribution. 



(a)Vt=Q (b)tft=156 (c)t/t=188 (d)t/t=219 (e}tft=453 




Figure 30: Time development of the dynamic Kelvin-Helmholtz instability for vr = 
0.98c > Vd- The density difference is scaled by n T = fi/g. The height in the lower 
figures represents the vorticity u> v . The numerical simulation was done under the periodic 
boundary condition along the x axis. The system size is 64£ x 64£ with £ = h/ \fgp. The 
time is scaled by r = ft/^.[Takcuchi, Suzuki, Kasamatsu, Saito and Tsubota: Phys. Rev. 
B 81 (2010) 094517, reproduced with permission. Copyright 2010 the American Physical 
Society] 

In the linear stage of the instability, a random seed, added in the initial 
stationary state, grows into a sinusoidal interface wave. The wave number 
of the sinusoidal wave corresponds to that of the unstable mode with the 
largest imaginary part |Im u\. As the amplitude of the wave becomes large, 
the sine wave is distorted [Fig. [30] (b)], and deforms into a sawtooth wave 
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[Fig. 1201 (c)]. Then, quantized vortices are released from the edges of the 
sawtooth waves. The vorticity u v is localized on the edges of the sawtooth 
waves and creates singular peaks [Fig. [30] (d)]. These peaks are released from 
the vortex sheet, becoming a singly quantized vortex with a circulation of 
K = h/m [Fig. [30] (e)]. The vorticity of the vortex sheet on the interface is 
reduced after the release of vortices, and then the velocity difference across 
the interface decreases locally. The relative velocity across the interface after 
the vortex nucleation is roughly estimated to be the total vorticity on the 
interface divided by the length of the sheet. Since six quantized vortices 
are released from the interface [Fig. [30] (f)], the relative velocity decreases 
by about 6k,/ L ~ 0.6c with the system size L = 64£ below the threshold 
V)> Then the instability stops, and vortex nucleation occurs no more. The 
released vortices continue to drift along the interface and the system does 
not recover the initial flat interface. 

Since the Kelvin-Helmholtz instability occurs locally around the inter- 
face, the relative motion of the two components does not ultimately vanish, 
in contrast to the countersuperflow instability demonstrated in the previous 
subsection. The instability can occur globally if the interface thickness be- 
comes comparable to the system size for a small inter-component interaction 
with gi2 ~ g. The instability phenomenon then becomes similar to that of 
countersuperflow instability. The crossover between counterflow instability 



and Kelvin-Helmholtz instability is investigated in Ref. [176 



4-4- 3. Nonlinear development of thermodynamic Kelvin-Helmholtz instabil- 
ity 

We next discuss the nonlinear time development of the Landau instability 
(thermodynamic Kelvin-Helmholtz instability) for vr > Vl- The dissipative 
dynamics can be qualitatively investigated by solving the dissipative GP 
equations, which are obtained by replacing the time-derivative term id t by 
(i - j)dt in the GP equations (ITUgj) Q- 

Figure I3T1 shows the time development of the instability obtained by solv- 
ing the dissipative GP equations. In the nonlinear stage, the interface has 
flattened troughs and peaked crests [Fig. l3~iT b)]. The patterns are diphyc- 
ercally asymmetric in contrast with the patterns of the dynamic Kelvin- 
Helmholtz instability. In this case, the vorticity is localized at the crests 
[Fig. l3TT c)] and a single quantized vortex is nucleated from each crest only 
into the upper side [Fig. IBTT d)]. 

The nucleated vortices are dragged away from the interface to the upper 
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direction due to the dissipation, which decreases the flow velocity of the 
2nd component (the blue region in Fig. 1ST]) at the rest fram e. The event 



of the vortex dragging is interpreted as phase slippage [1841 . The vortex 



nucleation stops after four vortices are nucleated, where the relative velocity 
across the interface decreases by Ak/L ~ 0.4c below the threshold Vl of the 
Landau instability. Then the interface recovers a flat form with less vorticity 
(relative velocity) than the initial state. 



(a) Vt^Q 



(b)tf r=844 (c)fr=912 (d) t/p=1250 (e) t/^1781 



4 J 





(n r n 2 )/n T 



-i o +1 



Figure 31: Time development of the thermodynamic Kelvin-Helmholtz instability for 
vr = 0.79c > Vl- We set the dissipation coefficient to be 7 — 0.03. 



4-4-4- Conclusion 

The interface modes are amplified due to dynamic instability when the rel- 
ative velocity across the interface exceeds a critical value in phase-separated 
two-component BECs. The instability is interpreted as the quantum coun- 
terpart of the Kelvin-Helmholtz instability when the interface motion is rep- 
resented using the hydrodynamic formalism. The nonlinear dynamics of the 
quantum Kelvin-Helmholtz instability are quite different from those of classi- 
cal Kelvin-Helmholtz instability governed by quantized vortices. The Landau 
instability for the interface modes in the presence of relative velocity, called 
the thermodynamic Kelvin-Helmholtz instability, causes vortex nucleation 
and phase slippage from the interface, which has no analog in classical fluid 
dynamics. The Kelvin-Helmholtz instability in trapped systems and the pos- 



sibility of its experimental realization was discussed in detail in Ref. [176 
We believe that the first observation of the dynamic Kelvin-Helmholtz insta- 
bility, the quantum counterpart of the classical Kelvin-Helmholtz instability, 
will soon be realized in future experiments. 
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5. Conclusions 

We have reviewed recent topics on quantum hydrodynamics (QHD) in 
superfluid helium and atomic BECs, chiefly focusing on the activity of our 
group. Quantized vortices were discovered in superfluid 4 He in the 1950s. 
However, they have recently grown in importance, for two reasons. The first 
reason is that the research of QT entered a new era since the mid 1990s 
leaving the previous studies almost limited to thermal counterflow. The 
second reason is the realization of atomic BECs in 1995. Modern optical 
techniques have enabled the direct visualization of quantized vortices, and 
multi-component BECs have further enriched the world of quantized vortices. 
In this concluding section, we describe the main motivation of this research 
and the interesting topics that are not addressed in the text. The discussions 
in this section are limited to the case at zero temperature where the normal 
fluid component is negligible. 

Comparing QT and CT reminds us of the motivation of studying QHD. 
Turbulence in a classical viscous fluid appears to be comprised of eddies. 
However, these eddies are unstable and not well defined. The circulation is 
not conserved and is not identical for each eddy. QT consists of a tangle of 
quantized vortices that have the same conserved circulation. Looking back 
at the history of science, reductionism, which tries to understand the na- 
ture of complex things by reducing them to the interactions of their parts, 
has played an extremely important role. The success of solid-state physics 
owes much to reductionism. In contrast, conventional fluid physics is not 
reducible to elements, and thus does not enjoy the benefits of reductionism. 
However, QT is different, being reduced to quantized vortices; reductionism 
is applicable to quantum turbulence. The main interests would be quan- 
tum hydrodynamic instability and QT beyond the instability. How can we 
approach these problems from reductionism! 

We should reveal the transition to QT and the nature of QT. Statistical 
quantities are useful in order to investigate the transition and the nature of 
QT. Here we list the possible statistical quantities. 

1. Energy spectra. The energy spectrum of fully developed QT is expected 
to obey the Kolmogorov law. However, the present understanding of 
the story is not so simple [7]. It is generally believed that there are two 
kinds of QT, namely quasi-classical turbulence and ultra-quantum tur- 
bulence. In quasi-classical turbulence, most of the turbulent energy is 
concentrated in the large-scaled eddies, typically on scales larger than 
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the mean inter-vortex distance I. This is quite similar to the case of 
classical homogeneous turbulence. The discreteness of each quantized 
vortex is not so relevant; they are expected to make some coherent 
structure (vortex bundles). Then the energy spectra would follow the 
Kolmogorov law. If we switch off the energy injection sustaining the 
turbulence, the vortex line density (VLD) decays as L oc t~ 3 / 2 . On 
the other hand, ultra-quantum turbulence has no quasi-classical mo- 
tion on scales greater than i. Most of energy is concentrated on scales 
smaller than t. When it decays, the VLD reduces as L oc t -1 . There 
are no direct observations of the energy spectra at such low tempera- 
tures. However, th e two typ es of decay L oc t~ 3 / 2 and t" 1 are observed 



experimentally [53l . Il85l Il86l |. The essential point is how to connect the 
configuration of quantized vortices with energy spectra. If a configura- 
tion of vortices is given, the energy spectrum is determined uniquely. 
Then, what types of vortex configuration can make quasi-classical or 
ultra-quantum turbulence? 
2. PDF(Probability density function) of superfluid velocity v s . This was 
discussed in Sec. 13.3.11 The PDF shows classical Gaussian in low 
velocity and non-classical power-law in high velocity. However, there 



are only a few theoretical works on this problem [45l . |46| . no system- 
atic studies in connection with energy spectra and the configuration of 
vortices. 

3. Vortex length distribution. If self-similar Richardson cascade appears, 
it is expected to yield some self-similar power-law in th e vo r tex l ength 



distribution. There are still a few numerical works |83|, I118L Ill9j . It is 
impossible to observe the distribution in superfluid helium, but possible 
in atomic BECs. 

Drag coefficient. This quantity was discussed in Sec. 13.41 The drag 
coefficient Co is inversely proportional to the velocity in laminar flow 
and of order unity in turbulent flow. This is a well-known story in CT, 



and confirmed in QT too [47j|. This change in can be another sign 
of transition to QT. 

It is possible and important to consider QT as a tra nsient state in the 



relaxation process far from thermal equilibrium |187 . 188] . QT corresponds 
to nonthermal fixed points in a nonperturbative quantum-field theoretic ap- 
proach, following some scaling law characteristic of dynamical critical phe- 
nomena. This kind of approach should make progress in near future. 
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Section H] described chiefly hydrodynamics and QT in two-comp onen t 



BECs. A spinor BEC is another important system for hydrodynamics [189 
Fujimoto and Tsubota investigated theoretically and numerically the GP 
model of spin-1 spinor BECs. They considered the two cases: one is the 
counterflow of two components with different magnetic quantum numbers in a 



uniform system 190] and the other is starting from a helical spin structure in a 
trapped system 1911 ] . When the interaction is ferromagnetic, the instability 
is amplified to spin turbulence in both cases, where the spectrum of the 
spin-dependent interaction energy exhibits a -7/3 power law, different from 
the Kolmogorov -5/3 law. This power law is understood from some scaling 
argument for the equation of motion of th e spin density vector. Since such 



spin density vector can be observed [192| . such spin turbulence could be 
realized and observed. 
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